Cédrick Lelionnais
email: cedrick.lelionnais@eseo.fr

Matthias Brun
email: matthias.brun@eseo.fr

Jérôme Delatour
email: jerome.delatour@eseo.fr

Olivier Henri Roux
email: olivier-h.roux@irccyn.ec-nantes.fr

Charlotte Seidner
email: charlotte.seidner@irccyn.ec-nantes.fr

Formal Behavioral Modeling of Real-Time Operating Systems

Keywords: Model Driven Engineering, Real-time operating systems, Behavioral modeling, Formalization, Time Petri Nets, Application deployment

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Nowadays, Real-Time Embedded Systems (RTESs) increasingly surround us in various domains (aircrafts, automotive, cell phones, robotics...). RTES engineers are confronted with the challenge of developing more complex, higher quality systems, with shorter development cycles at lower costs. Within this context, reuse, maintainability and portability have become major issues in RTES design processes.

From this point of view, some frameworks have been developed with the aim of generating code. Starting from a detailed model written in such a language as AADL (Society of Automotive Engineer (SAE), 2004) or UML (Object Management Group (OMG), 2007b), the specific code for the deployment of an application on a Real-Time Operating System (RTOS) is then generated. These code generating frameworks are generally used after some Verification and Validation (V&V) activities have been performed on the models. In this way, some research works have dealt with the transformation of those models into formal models, using such formalisms as Time Petri Nets (TPNs) or Finite State Machines (FSMs)...) as presented for instance in [START_REF] Renault | From aadl architectural models to petri nets: Checking model viability[END_REF] [START_REF] Berthomieu | Formal Verification of AADL models with Fiacre and Tina[END_REF].

However, the correctness of the application deployed on the RTOS is not verified during these V&V activities. It is therefore difficult to ensure that the generated code effectively corresponds to what was designed. Moreover the consideration of the RTOSs is often embeded within the processus of generation (or seen in simplified terms). This approach is then far from flexible since each translation must be adapted in accordance with the targeted RTOS. Due to the significant number of existing RTOSs, this results in having to maintain a large amount of translations.

In consequence, this paper relates to the following issues. Which methodology should be adopted in order to both deploy an application on a targeted RTOS and check the correctness of such a deployment? Is it possible to consider any RTOS during the deployment? (and if so, how?) How can the behavior of any RTOS be described within such frameworks? To bring an answer to these questions, a method has been experimented to take these factors into account within a single framework. The basic idea is to add the formal description of any RTOS. Each description can then be considered as a parameter of the transformation, so as to keep the framework generic. Lastly, this formalization should let V&V tools able to verify the correctness of such deployments.

This paper is divided into the following sections. Section 2 presents some works in conjunction with the issues described above. As a result of that research, the approach adopted is developed in Section 3. It has proven the need to consider the behavior of RTOSs during a deployment of an application. Then Section 4 presents the strategy proposed to formally take into account the behavioral description of RTOSs. Consequently the benefits and the limits of the strategy are showed in Section 5 and discussed in Section 6. Finally we conclude in Section 7.

RELATED WORK

Code generating frameworks

Code generating frameworks are used in a context of application deployment. Code generation is the process of transforming code from one representation to another one. Often, this is from a higher level, such as an UML diagram modeling an application, to a lower level, such as a C code program that is portable to an execution software platform.

From a real-time point of view, code generating frameworks are used in a context of application deployment to a specific RTOS. The information concerning RTOSs are given in a program which implements the deployment. For instance, Ocarina described the information of a middleware with Ada code in its Gaia generator [START_REF] Vergnaud | Ocarina, a compiler for the AADL[END_REF].

Similarly, the Real-Time Workshop (RTW) tool (The MathWorks, 2007) uses templates catching both technological concepts of the RTOS and programming language information. These templates structure the source code to generate.

As a result, these frameworks involve an implicit description of the RTOSs during the processus of generation. This raises the problem of genericity of these frameworks since each deployment corresponds to a unique code generator.

RTOSs modeling

In order to describe the Operating Sytems (OSs) concepts, some model-based frameworks have been developed. Some works have focused on both structural and behavioral aspects of the OSs. Revolving around these two axes, the description of the OSs is made explicit and observable through formal models. Explicit means that such a description is entirely dedicated to the representation of the OS.

Metropolis [START_REF] Team | The metropolis meta model -version 0.4[END_REF] uses formal languages such as the Linear Temporal Logic. It offers the possibility to verify some properties of an application deployed on a wide range of OSs.

Generic Modeling Environment (GME) [START_REF] Davis | GME: the Generic Modeling Environment[END_REF] contributes to the description of OSs which are not necessarely real-time and analysis platforms. GME proposes assistance tools for the migration from one OS to another. In addition, this environment offers the possibility to insert formal metamodels. The descriptions thereby made can be used by V&V tools.

Nevertheless, the genericity of the modeling languages used within these frameworks does not make the modeling of platforms easier dedicated to a particular domain, such as the real-time domain, any easier.

More specifically to the RTESs, Ptolemy [START_REF] Lee | Overview of the Ptolemy project[END_REF] is a framework which describes execution models. Most of these models support the actor-oriented design focused on the concurrency and the communication between the components of a RTES. However, the concepts of structural representation from this approach are only intended for the applications modeling, but not for RTESs themselves.

In spite of the fact that the OSs description is explicit, the works outlined above are less in line with the RTESs. Furthermore, the separation of the concerns does not clearly appear. This impacts on the intervention of each domain specialist such as the application or the OSs technology. The lack of genericity previously noted arises once again.

Other researches have led to consider the description of RTOS during the application deployment. This approach relies on a Model Driven Engineering (MDE) context. MDE promotes the distinction of platform-specific modeling artefacts from those that are platform-independent. The platform (RTOS) model is considered as a parameter of this deployment. Each model is described by a modeling language. Concerning the RTOS modeling, we can note the Software Resource Modeling (SRM) [START_REF] Thomas | Embedded Systems Specification and Design Languages[END_REF] package which is a MARTE's UML profile (Object Management Group (OMG), 2007a). With SRM, RTOSs can be modeled using stereotyped concepts from the real-time software domain.

For another example, RTEPML (Real-Time Embedded Platform Modeling Language) [START_REF] Brun | Contribution on the software execution platform integration during an application deployment process[END_REF] was developped with the aim of defining concepts dedicated to the real-time domain for modeling RTOSs.

Contrary to the first noted frameworks, SRM and RTEPML better separate the concerns previously highlighted (RTOSs structure, transformations, deployments choice...). Nevertheless, only the structural aspect is taken into account with these modeling languages. The behavior resulting from the concurrent activities of the RTESs components is not observable in the deployed application models. As a result, no V&V activities can be applied to check the correctness of such a deployment.

APPROACH ADOPTED

The previous section showed the interest to well distinguish the concerns for the intervention of each specialist for real-time applications deployment. In addition, it has been raised to explicitly consider the description of the RTOSs behavior within a code generating framework. Our approach relies therefore on the extension of a model-based framework to formally observe the behavior of a deployed application on a RTOS. This should add other concerns such as the RTOSs behavior and some V&V activities on the deployment at design phase. As a result, we have chosen the RTEPML framework as it has been developed within our team. In addition this modeling language already contributes to the explicit description of the RTOSs. However, the behavioral aspect of RTOSs can not be modeled with RTEPML. Therefore modeling behavior with RTEPML means extending this language. That is why we have sought to enrich it in order to define behavioral concepts of RTOSs.

RTEPML with MDE

Based on the MDE initiative, the application deployment emanates from semi-generic transformations of models. MDE is based on an architecture (Model Driven Architecture (MDA)) (Object Management Group (OMG), 2001) (see Figure 1) involving a design process where a Platform-Independent Model (PIM) (description of the application without considering the platform) is transformed into a Platform-Specific Model (PSM) (description of the application deployed on the platform), according to a Platform Description Model (PDM) (description of the platform). Each model conforms to a modeling language, also named metamodel.

In accordance with Figure 1, RTEPML is a metamodel which has been built to model PDMs. The concepts of the targeted RTOS are then defined by a PDM. So far, the modeling with RTEPML goes through the structural characterizations (resources Figure 1: MDA principle and services) of RTOSs. The resources represent the concepts offered by the platform. For instance, in a RTOS, the tasks, the semaphores, the mailboxes are considered as resources. The services are offered by the resources and caught at API's level (Application Programming Interface). For example, services associated with semaphores are the take and release primitives.

Consequently, in conjunction with this approach, the deployment consists both in instantiating RTOS resources and in using RTOS call services for implementing the application. Then the generated PSM is provided for a code generation process. Finally, it is worth noting that the transformation serving to the deployment is generic thanks to the explicit consideration of the RTOSs.

Illustration around semaphore sharing

We have mentioned the idea to extend RTEPML at the beginning of this section. This extension should enable to consider the RTOS behavior during an application deployment. To illustrate this behavioral lack, a deployment is presented Figure 3 1 based on an example of Robot application (see Figure 2). Obviously several concepts instantiations of the considered RTOS are possible following the application. Accordingly, this illustration is focused on a particular case : the semaphore sharing.

In this way, the PIM representing the robot application combines two periodic activities (a Driver activity for driving the robot and a Vigil activity for avoiding barriers). On the one hand, the Driver in- Then two wheels allow the robot to move thanks to two settings updated according to the activities. Finally each periodic activity owns a priority according to its role.

The PDM described by RTEPML for this deployment is a model of the nxtOSEK's RTOS [START_REF] Nxtosek | NXTOSEK Operating System[END_REF]. nxtOSEK is based on OSEK/VDX standard (OSEK/VDX [START_REF] Group | OSEK/VDX Operating System Specification[END_REF] and is widely used in LEGO MINDSTORMS NXT educational projects. According to our needs, some resources have been defined to be instantiated with regards to the application. The Task concept is instantiated to represent each activity, the Function concept to execute each program and the Resource (semaphore concept in nxtOSEK) concept to protect each wheel setting updating.

The PSM results in a model of the application deployed on the nxtOSEK platform. However the concurrent activities which are governed by both application and technological concepts do not concretely appear in the final model. For instance we do not know whether each Resource (protection for updating each wheel setting) is available or not (empty or full). Moreover the state in which each Task (the Driver T and the Vigil T) is set cannot be clearly distinguished. Finally the impact of the properties such as the priority and the period of the tasks is not directly verifiable.

Consequently this illustration shows how difficult it is to interpret the behavior of the resources instantiated. The behavioral aspect of the generated model is not observable. In this approach, it is therefore difficult (if not impossible) to undertake some V&V activities to verify the correctness of the deployment. In order to be able to check the properties on the deployed application, formal behavioral aspects have to be introduced within the description of the software platform.

STRATEGY OF BEHAVIORAL FORMALIZATION

In accordance to the previous conclusion, RTEPML has been enriched to bring a behavioral aspect in the RTOSs description. On this point, we have extended the abstract syntax of RTEPML in which executive concepts have been defined to explicitly describe any platform.

Once RTEPML extended (newly called RTEPML BEHAVIOR), the implementation of a transformation process could have been put in place to make the model of the deployed application behavioral. To meet the requirements of V&V activities, the formalization of the PSM is necessary. RTEPML BEHAVIOR has therefore been designed to take into account any formal language. This solution has been thought in order to make the transformation independant of other formal languages.

Behavioral Consideration

Taking into account the RTOSs behavior necessitates the identification of the behavior of each resource (scheduling task, state of a semaphore...). The same applies to the behavioral description of each service offered by resources (task terminating, semaphore acquisition...). Lastly, each behavioral impact associated with each resource property (task priority, task period, semaphore tokens number...) should also be treated during the deployment.

To keep on the previous focus, an example of behavioral consideration around the semaphore is depicted Figure 4. As explained earlier, some concepts have been added in RTEPML BEHAVIOR to represent the behavior of the RTOSs resources. Thus, the behavioral prototype of the semaphore concept, called ResourceBehavioralPrototype, has been defined on the Behavioral PDM. As a reminder, Resource (on Behavioral PDM) represents the semaphore concept of nxtOSEK RTOS, whereas "Resource" is the RTEPML BEHAVIOR metaconcept allowing to describe any RTOS resource. In other words, Resource conforms to MutualExclusionResource.

The formalism used to express the behavioral prototype is Time Petri Net (TPN) (see [START_REF] Boyer | On the compared expressiveness of arc, place and transition time Petri nets[END_REF] for a survey). TPNs have been prefered because of their expressiveness. This formalism is indeed adapted to the description of concurrency activities such as resources schedulability, resource sharing, call services sequence, time constraints... Furthermore, TPNs propose known time extensions (such as stopwatch TPNs) which allow to model real-time Figure 3: Robot application deployment schedulings implemented in V&V tools like Romeo [START_REF] Gardey | Roméo: A tool for analyzing time Petri nets[END_REF] [START_REF] Lime | Romeo: A parametric model-checker for Petri nets with stopwatches[END_REF] or Tina [START_REF] Berthomieu | The tool tina -construction of abstract state spaces for Petri nets and time Petri nets[END_REF]. As can be seen on Figure 4, a behavioral prototype can be composed of several behavior elements. Indeed, a RTOS resource includes services and properties. The appearance of roles (defined at metamodel level) within the Behavioral PDM, enables to precise certain behavior elements. In the nxtOSEK semaphore case, serviceBehavioralDescription role comes to describe the services behavior of getResource and releaseResource, whereas defines role comes to define the value of the semaphore.

Deployment with RTEPML BEHAVIOR

The behavioral description of RTOSs has been made possible with RTEPML BEHAVIOR. Moreover the possibility to define formally the behavior elements, contributes to the formalisation of the application deployment. This important point should enable to verify its correctness before the code generation.

In accordance with the adopted approach, a transformation has been developed to consider explicitly a behavioral PDM of the targeted RTOS. We have taken the same robot application from the previous section to illustrate an example of formalisation of the deployment (see Figure 5). For the sake of simplicity and clarity, we have reduced the Robot PIM to the Driver T instance representing the periodic Driver activity. Driver T has an entry point on the Moving Strategy in which only a service is called for terminating itself.

The involved prototypes in the transformation are those of the resources for the needs of the application. Once the prototypes are detected, the latter are cloned according to each instance of the application. It may be observed that the cloned prototypes corresponding to the behavior of the Driver T and the Moving Strategy (including the service of terminating) have been composed through this example.

As a result, the TPN Behavioral PSM achieved is more expressive than the previous one. Its expressiveness helps us to easily see its liveness (the different states of Driver T, its periodicity...).

The composition of the building blocks represented by cloned behavioral prototypes (semaphore, task, communication, priority...) make up the generated model. In TPNs case, the composition is carried out by merging places interfacing all the blocks. This presents some advantages. The semantics of the places merging is obvious, while the one of the transitions merging is ambiguous. Indeed, the fact of assigning different time constraints on mergeable transitions makes their merging difficult. In addition, the In this paper, the purpose is neither to describe exhaustively the modeling of the behavioral prototypes and the composition of their clones, nor to prove that it is correct. The aim is only to present the approach and the methodology which allows it. That is why no proof of the building validity has been given. Notwithstanding, this proof is greatly simplified by our composition approach.

BENEFITS OF THE STRATEGY

The main benefits of our methodology are both the formalization and the moving of the application deployment into the design phase. This should allow to ensure its correctness before generating code. Basing our approach on MDE, a formalization of the deployed application model can be carried out. V&V activities will become applicable thanks to the use of model-checking tools.

A second benefit is the separation of the skills highlighted in the second section. The distinction of modeling artefacts allows both RTOSs behavior and formal verification experts to interact on the application deployment (see Figure 6). Thanks to RTEPML BEHAVIOR, the specialist of the RTOSs behavior can be assisted by the one more highly skilled in the formal verification. The latter one can also interpret the most adapted formalism following the case.

To put our strategy into practice, we have chosen to verify a necessary and well-known condition in the Despite the combinatorial explosion risk of the generated TPN, ROMEO has given a deadlock outcome. This diagnosis could seem foreseeable from the specification since the wheels settings updating is protected by semaphores. Notwithstanding, this verification helps us to evaluate the right protocol to use for this application. Indeed, the Priority Ceiling Protocol (PCP) (OSEK/VDX [START_REF] Group | OSEK/VDX Operating System Specification[END_REF] should have been chosen to remedy this situation during the deployment.

We have exposed the deadlock case but other behavioral properties could have been tested. For instance we could deal with liveness properties such as termination of actions, occurency of expected events... But we could also handle time constraints such as Worst-Case Execution Time (WCET), respect of a deadline...

LIMITS

This experimentation has showed the feasability of generating a formal model to observe the behavior of the application deployed on a RTOS. Nonetheless, the behavior of some of concepts (the notification resources such as the events or the communication resources such as the messages...) are not yet taken into account within RTEPML BEHAVIOR. As a result, an optimization of the latter modeling language must be widened in order to cover such concepts.

Another limit raised is the interpretation of other formalisms to model the behavior during the deployment. Only TPNs have been presented through our strategy but the implementation of other formal languages deserves to be examined.

Lastly, a deadlock has been detected during the deployment. It has been deduced that the deployment was not correct. That being so, this decision does not affect the correspondence between the application and the targeted RTOS. Indeed, a protocol changement (Priority Ceiling Protocol (nxtOSEK, 2009) (OSEK/VDX [START_REF] Group | OSEK/VDX Operating System Specification[END_REF]) is just sufficient to prove it. This leads us to say that the framework does should integrate decision support tools for guiding the experts involved in the RTOS choice.

CONCLUSION

Through this study, a methodology has been presented to consider the behavior of a real-time execution software platform (RTOS) during an application deployment. The deployment which enables to implement an application on a specific RTOS, has been moved at design phase. The purpose of this choice is to enable the different specialists to intervern more easily on the deployment following its domain. Indeed, thanks to the modeling with RTEPML and the MDE approach, the separation of concerns (application and RTOS) has been made explicit.

The description of the RTOSs behavior has been made feasible by enriching RTEPML (RTEPML BEHAVIOR). In addition, each behavioral model can be formally described.

A transformation has been developed to generate a formal model of the deployed application. The formalization has the advantage of applying V&V activities for checking the correctness of the deployment, even before generating the code useful to the implementation of the application.

However, certain limits of our methodology have arisen. Future prospects are the subject of our next works to improve and to check the feasibility of such a strategy. In accordance to the previous section, the behavioral concepts missing in RTEPML BEHAVIOR will have to be created. Then, to ensure the genericity of the formalisation, other languages will have to be interpreted. Finally, a large number of formal proofs will have to be written within our framework. This should contribute to the verification of necessary and sufficient conditions for validating such deployments.

Figure 2 :

 2 Figure 2: Model of a Robot application

Figure

  Figure 4: MER behavior

Figure 5 :

 5 Figure 5: Formalisation of robot application in TPN

Figure 6 :

 6 Figure 6: Intervention of specialists

The languages used in Figures

to 5 come from the graphic description of each metamodel (their concrete syntax)

http://romeo.rts-software.org


Jean-Marc Robert

Software Implementation of Parallelized ECSM over Binary and Prime Fields

Keywords: Elliptic curve cryptography, parallel algorithm, efficient software implementation

Recent developments of multicore architectures over various platforms (desktop computers and servers as well as embedded systems) challenge the classical approaches of sequential computation algorithms, in particular elliptic curve cryptography protocols. In this work, we deploy different parallel software implementations of elliptic curve scalar multiplication of point, in order to improve the performances in comparison with the sequential counter parts, taking into account the multithreading synchronization, scalar recoding and memory management issues. Two thread and four thread algorithms are tested on various curves over prime and binary fields, they provide improvement ratio of around 15% in comparison with their sequential counterparts.

Introduction

Elliptic curve cryptography (ECC) is widely used in a large number of protocols: secret key exchanges, asymmetric encryption-decryption, digital signatures... The main operation in these protocols is the scalar multiplication (ECSM) defined as k • P where P is a point of order r on an elliptic curve E(F q ) and k ∈ [0, r[ is an integer. The scalar multiplication is computed with Double-andadd approaches which consist of sequences of several hundreds of doublings and additions of curve points. It is thus a costly operation which might be implemented efficiently.

In this paper we consider parallel approaches for software implementation of scalar multiplication. There are two versions of the Double-and-add scalar multiplication: the left-to-right and the right-to-left depending on the way the bits of k are scanned. On the one hand, the left-to-right version cannot be parallelized due to the strong dependence of the consecutive doublings and additions. On the other hand, the right-to-left version is easier to parallelize: this was noticed by Moreno and Hasan in [START_REF] Moreno | SPA-Resistant Binary Exponentiation with Optimal Execution Time[END_REF]. Indeed, in [START_REF] Moreno | SPA-Resistant Binary Exponentiation with Optimal Execution Time[END_REF], the authors provide an algorithm consisting in one thread producing the points 2 i P through consecutive doublings, which are then consumed by a second thread performing all the necessary additions. They did not provide any implementation results of their approach. In practice this can be challenging to implement efficiently the synchronizations between the two threads.

When the elliptic curve is defined over a binary field F 2 m , a formula exists (cf. [START_REF] Knudsen | Elliptic Scalar Multiplication Using Point Halving[END_REF][START_REF] Fong | Field Inversion and Point Halving Revisited[END_REF]) which computes efficiently the halving of a point, i.e., 1 2 P . This makes possible to perform the scalar multiplication through a sequence of halvings and additions of points. This can be used to parallelize the scalar multiplication into two totally independent threads: one thread performing a halve-and-add scalar multiplication and a second thread performing a double-and-add. This approach has been implemented by Taverne et al. in [START_REF] Taverne | Speeding Scalar Multiplication over Binary Elliptic Curves using the New Carry-Less Multiplication Instruction[END_REF] showing a significant speed-up compared to non-parallelized versions.

In this paper we first explore the implementation of the two threads parallel approach of Moreno and Hasan [START_REF] Moreno | SPA-Resistant Binary Exponentiation with Optimal Execution Time[END_REF]. Specifically, we analyze three different strategies to perform synchronization between both threads: using signals, mutexes or busy-waiting approaches, we propose a synchronization strategy based on this analysis. We also study the best approach for the coding of the integer k: this impacts the number of additions and post-computations, i.e., the work load of the thread performing the additions.

We then investigate a four thread parallelization of the scalar multiplication in E(F 2 m ). This approach combines the Double/halve-and-add algorithm of [START_REF] Taverne | Speeding Scalar Multiplication over Binary Elliptic Curves using the New Carry-Less Multiplication Instruction[END_REF] with the approach of Moreno and Hasan.

We provide experimental results for two curves defined over a prime field p = 2 255 -19 and for the two binary elliptic curves B409 and B233 recommended by NIST in [START_REF] Director | Digital Signature Standard (DSS)[END_REF]. Our experimental results show that the parallelized scalar multiplication is up to 15 % faster than their non-parallelized counterparts (depending of the curve type and the field size).

The remaining of the paper is organized as follows: in Section 2 we review basic definitions of elliptic curve and scalar multiplication algorithms. In Section 3, we present our implementation approaches of scalar multiplication. We then provide in Section 4 the experimental results and comparisons with the state of the art. We end the paper in Section 5 with some concluding remarks.

Background on elliptic curve scalar multiplication

In this section, we briefly review basic results concerning elliptic curve and their use in cryptography. For further details on this matter we refer the reader to [START_REF] Hankerson | Guide to Elliptic Curve Cryptography[END_REF]. An elliptic curve over a finite field E(F q ) is the set of point (x, y) ∈ F 2 q satisfying a smooth curve equation of degree 3 in x and y plus a point at infinity O. A group law can be defined using the so-called chord-and-tangent approach, providing formulas in terms of point coordinates which compute doubling 2P and addition P + Q in the group. The element O is the neutral element of the group. Cryptographic protocols are based on the intractability of the discrete logarithm problem: given a generator of the group P and a point Q, compute k such that Q = kP . The most costly operation involved in most ECC protocols is the scalar multiplication: given P ∈ E(F q ) and an integer k, the scalar multiplication consists in computing kP = P + P + • • • + P (k times). The elliptic curves used in practice are defined either over prime field F p with p prime or over binary field F 2 m . In the remainder of this section, we briefly review explicit formulas and algorithms for scalar multiplication over these two fields.

Scalar multiplication over prime field

Weierstrass Elliptic curve. An elliptic curve E over a prime field F p is generally defined by a short Weierstrass equation:

E : y 2 = x 3 + ax + b, (a, b) ∈ F 2 p .
Then, in this case, addition and doubling on E(F p ) works as follows: let P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), P 3 = (x 3 , y 3 ), be three points of E such that P 3 = P 1 +P 2 , then we have:

x 3 = λ 2 -x 1 -x 2 , y 3 = λ(x 1 -x 3 ) -y 1 , where λ = y2-y1 x2-x1 if P 1 = P 2 , λ = 3x 2 1 +a 2y1 + x 1 if P 1 = P 2 .
Jacobi Quartic curves over prime field. This curve was suggested by Billet et al. in [START_REF] Billet | The Jacobi Model of an Elliptic Curve and Side-Channel Analysis[END_REF]. The curve equation of E is:

y 2 = x 4 - 3 2 θx 2 + 1, θ ∈ F p .
For such curve, the addition and doubling formulas are unified. Let P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), P 3 = (x 3 , y 3 ), be three points of E such that P 3 = P 1 + P 2 , then we have:

x 3 = (x 1 y 2 + y 1 x 2 )/(1 -(x 1 x 2 ) 2 ), y 3 = ((1 + (x 1 x 2 ) 2 )(y 1 y 2 + 2ax 1 x 2 ) + 2x 1 x 2 (x 2 1 + x 2 2 ))/(1 -(x 1 x 2 ) 2 ) 2 .
The Jacobi Quartic curve is isomorphic to the following Weierstrass elliptic curve: y 2 = x 3 + ax + b where a = (-16 -3θ 2 )/4 and b = -θ 3 -aθ.

Elliptic curve point operations. The most expensive field operation is the inversion which roughly requires several tens of field multiplications. In order to avoid such operation, additions and doublings utilize projective coordinate system. In our implementation, we consider two systems: the Jacobian coordinate where the point (X : Y : Z) corresponds to the affine point (X/Z 2 , Y /Z 3 ) and the XXY ZZ coordinate system where the point (X : XX : Y : Z : ZZ) corresponds to the affine point (X/Z, Y /ZZ) with XX = X 2 and ZZ = Z 2 . Explicit formulas for addition and doubling in these systems can be found in [START_REF]Explicit formula database[END_REF] The resulting complexities are shown in Table 1, which shows that the complexities of the Jacobi Quartic curve operations are better than for the Weierstrass equation case. Moreover, based on the elliptic curve formula database in [START_REF]Explicit formula database[END_REF], the Jacobi Quartic curves provide the most efficient point operation among all known curves and formulas. This is the reason why we used such curve and these formulas in our implementations. 1. Weierstrass curve and Jacobi Quartic curve point operations, M = multiplications, S = squaring, R = field reduction.

Scalar multiplication algorithm. The basic method to compute a scalar multiplication consists in scanning the bits

k i of k = t-1 i=0 k i • 2 i
and performing a sequence of doubling followed by an addition when k i = 1. This approach is described in Algorithm 1.

In order to reduce the number of additions, the non adjacent form (NAF) and the window non adjacent form (W-NAF) recoding of the scalar are wellknown methods, which reduce the number of non zero digit representing the scalar. In the binary scalar representation, half of the digits are either zero or one on average. In the NAF representation, one uses three digits instead of two:

k = t i=0 k i • 2 i with k i ∈ {-1,
0, 1} and there are only t/3 non zero digits k i on average.

The W-NAF representation extends this concept by using more digits:

k = t i=0 k i • 2 i with k i ∈ {-(2 w-1 -1), .
. . , -5, -3 -1, 0, 1, 3, 5, . . . , (2 w-1 -1)}. The number of non zero digits is now t/(w + 1) on average. Algorithm 1 can be adapted to use k recoded as NAF or W-NAF. The complexities of the resulting scalar multiplication are given in Table 2.

nb. of doublings nb. of additions Double-and-add t -1 t/2 NAF Double-and-add t t/3 W-NAF Double-and-add t t/(w + 1) + 2 w-2 -1 Table 2. Complexity comparison between binary, NAF and W-NAF scalar representation in terms of t the bit length of the scalar.

The reader may refer to [START_REF] Hankerson | Software Implementation of Elliptic Curve Cryptography over Binary Fields[END_REF] for further details and algorithms to compute NAF and W-NAF representation.

Elliptic curve scalar multiplication over binary field

An elliptic curve E over a binary field F 2 m is the set of points P = (x, y) ∈ F 2 2 m satisfying the following equation:

E : y 2 + xy = x 3 + ax 2 + b, (a, b) ∈ F 2 2 m .

Algorithm 1

Left-to-Right Double-and-add

Require: k = (kt-1, . . . , k1, k0), P ∈ E(F 2 m ) Ensure: Q = k • P 1: Q ← O 2: for i from t -1 downto 0 3: Q ← 2 • Q 4: if ki = 1 then 5: Q ← Q + P 6: endif 7: endfor 8: return (Q) Algorithm 2
Right-to-left Halve-and-add

Require: P ∈ E(F 2 m ) of order r and k an integer in [0, r[ Ensure: Q = k • P 1: Compute k = 2 t • k mod r = t i=0 k i 2 i with t = log 2 (r) + 1 2: Q ← O 3: for i from t downto 0 4: if ki = 1 then 5: Q ← Q + P 6: endif 7: P ← P/2 8: endfor 9: return (Q)
Let P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), P 3 = (x 3 , y 3 ), be three points of E such that P 3 = P 1 + P 2 , then we have:

x 3 = λ 2 + λ + x 1 + x 2 + a, y 3 = (x 1 + x 3 )λ + x 3 + y 1 ,
where

   λ = y1+y2 x1+x2 if P 1 = P 2 , λ = y1 x1 + x 1 if P 1 = P 2 .
(

) 1
Elliptic Curve Scalar Multiplication with halving. It was noticed by Knudsen in [START_REF] Knudsen | Elliptic Scalar Multiplication Using Point Halving[END_REF] that over a binary field, halving of points is possible in case of points of odd order since 2 admits an inverse modulo the order of the point. In other words, point halving is the reciprocal operation of point doubling: given

Q = (u, v) ∈ E(F 2 m
), one looks for P = (x, y) ∈ E(F 2 m ), P = -P such as Q = 2 • P . Based on equation (1), we know that x, y, u and v satisfy the following relations:

λ = x + y/x (2)
u = λ 2 + λ + a (3) v = x 2 + u(λ + 1) (4)
Consequently, in order to compute P , we first have to solve equation (3) to get λ (which means solve λ 2 +λ = u+a), then, equation ( 4) gives x = v + u(λ + 1), and finally, equation ( 2) gives y = λx + x 2 . The reader may refer to Knudsen in [START_REF] Knudsen | Elliptic Scalar Multiplication Using Point Halving[END_REF] and Fong et al. in [START_REF] Fong | Field Inversion and Point Halving Revisited[END_REF] for further details. In practice, this can be implemented efficiently and has roughly the same cost as two field multiplications (see [START_REF] Taverne | Speeding Scalar Multiplication over Binary Elliptic Curves using the New Carry-Less Multiplication Instruction[END_REF]).

The Double-and-add method can be modified into a Halve-and-add scalar multiplication. Preliminary, we need to change the scalar. Assuming the point P to be multiplied is of odd order r,

we compute k = 2 t • k mod r = t i=0 k i 2 i with t = log 2 (r) + 1. Then, we have k ≡ k /2 t ≡ t i=0 k i 2 i-t mod
r and the scalar multiplication can be computed as follows:

k • P = (k t + k t-1 • 2 -1 + . . . + k 0 2 -t ) • P.
This can be computed as a sequence of halvings and additions as shown in Algorithm 2.

Cost of elliptic curve point operations. Over a field of characteristic 2, and in order to avoid the inversions during the computation, which is the most expensive field operation again, one may use projective coordinate systems. The most interesting systems are the Lopez-Dahab (LD, as shown in [START_REF] Hankerson | Software Implementation of Elliptic Curve Cryptography over Binary Fields[END_REF]) and the Kim-Kim (KK, see [START_REF] Kim | A New Method for Speeding Up Arithmetic on Elliptic Curves over Binary Fields[END_REF]) projective coordinate systems. With such point representation, the addition and doubling operations do not include any inversion as shown in Table 3, and the whole scalar multiplication is computed with a significant speed-up. Table 3 shows that the complexities of KK are slightly better and then, when possible, we give the preference to the KK coordinate system. 3. Elliptic curve point operations, M = multiplications, S = squaring, SR = square root, QS = quadratic solver, R = reduction.

Strategies for parallel implementation of scalar multiplication

In this section, after a quick review of the implementation strategies used for the field operations, we expose how we elaborate the parallelized algorithm, taking into account all the constraints for such concurrent programming.

The platform used for the experimentations is an Optiplex 990 DELL , with a Linux 12.04 operating system. The processor is an Intel core i7 -2600 Sandy Bridge 3.4GHz. This processor owns four physical cores, which corresponds to the maximum thread number of our implementations. The code is written in C language and compiled with gcc version 4.6.3.

Field implementation strategies

Prime field implementation strategies. We considered the prime field F p , with p = 2 255 -19, which was introduced by Bernstein in [START_REF] Bernstein | Curve25519: New Diffie-Hellman Speed Records[END_REF]. To compute the field operations, we reused the publicly available code of Adam Langley in [START_REF] Langley | C25519 code[END_REF]. Based on our experiments, the code of Langley is significantly more efficient compared to low level functions of the GMP library [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF] for the considered field. In the code of Adam Langley a field element is stored in a table of five 64 bit words, each word containing only 51 bits. This allows a better management of carries in field addition and subtraction operations. The field multiplications and squarings are performed in two steps, which are multiprecision integer multiplication (respectively squaring) and a modular reduction of integer of 510 bit size (less than p 2 ) into field element of size 255 bits (reduction modulo p). The multiprecision integer multiplications and squarings are computed with the schoolbook method. The squaring operation is optimized with the usual trick which reduces the number of word multiplications. The reduction modulo p = 2 255 -19 of 510 bit size integer consists in multiplying by 19 the 255 most significant bits and adding the result to the 255 least significant bits. An inversion of a field element is computed using the Itoh-Tsujii method [START_REF] Itoh | A Fast Algorithm for Computing Multiplicative Inverses in GF(2 m ) Using Normal Bases[END_REF]: a -1 ≡ a p-2 mod p, and the exponentiation to p -2 is performed with a sequence of squarings and multiplications.

Binary field implementation strategies. Our implementations deal with NIST recommended fields F 2 233 = F[x]/(x 233 + x 74 + 1) and F 2 409 = F[x]/(x 409 + x 87 +1). Concerning the binary polynomial multiplication, we apply a small number of recursions of the Karatsuba algorithm. The Karatsuba algorithm breaks the m bit polynomial multiplication into several 64 bit polynomial multiplications. Such 64 bit multiplications are computed with the PCLMUL instruction, available on Intel Core i7 processors. Due to the special form of the irreducible polynomials, the reduction is done with a small number of shifts and bit-wise XORs. We compute the field inversion with the Itoh-Tsujii algorithm, that is a sequence of field multiplications and multisquarings performed with look-up table. For field squaring, square root and quadratic solver (needed in halvings), we also use a look-up table method, which is the fastest way according to our tests.

Remark 1. The use of Karatsuba for binary fields and schoolbook method for the prime field is due to the relative cost of word addition compared to word multiplication and to carry managements on our platform. Indeed, integer word additions and multiplications have roughly the same cost (1 vs 2 cycles). The use of Karatsuba algorithm for F p decreases the number of word multiplications, but, in counter part, it increases the number of additions and carry managements. For binary field the relative cost of addition (bit-wise XOR) and multiplications (PCLMUL instruction) is more important: 1 cycle vs 10 cycles. In this case Karatsuba is efficient to decrease the timing of a field multiplication.

Parallelization

The left-to-right Double-and-add algorithm (see Algorithm 1 page 5) does not allow any parallelization of the computations, due to the read-after-write dependency inside each loop iteration, between step 5 (addition) and step 3 (doubling). It is necessary to use the right-to-left variant of this algorithm (see Algorithm 3) which allows the parallelization. Indeed Algorithm 3 can be parallelized into two threads as follows:

-A producer-thread performing the sequence of doublings generating the points 2 i P . -An addition-thread accumulating the points generated by the producerthread.

In the sequential case, the left-to-right Double-and-add algorithm (Algorithm 1) is better, because the point addition in step 5 can use a mixed coordinate addition. This is faster than the projective addition used in the right-to-left version in step 4 (Algorithm 3). We will see that this penalty is overcome in most of the cases, thanks to the parallelization. The Halve-and-add algorithm (Algorithm 2 page 5) can also be parallelized with two threads. Indeed, since the computation in step 7 of Algorithm 2 only depends on the same step in the previous loop iteration (read-after-write dependency), the sequence of halvings (step 7) can be performed in a separate thread (the producer-thread) and the addition in an addition-thread which accumulates the points generated by the producer-thread.

Algorithm 3 Right-to-left Double-and-add

Require: k = (kt-1, . . . , k1, k0), P ∈ E(Fq) Ensure: The safest way to guarantee absolute correct computation is to use a strong synchronization device, processing the computation by small batches: the producerthread computes and stores a small batch of point doublings/halvings, sends a signal in order to trigger the addition computation in the addition-thread only concerning the batch in shared memory. In parallel, the producer-thread goes on with the next batch and the addition-thread waits the end of each batch before processing the corresponding additions (in the way described by Mueller in [START_REF] Mueller | A Library Implementation of POSIX Threads under UNIX[END_REF] or by Tannenbaum in [START_REF] Tannenbaum | Modern Operating Systems[END_REF]). In our case, on the one hand, the batch size has to be small to compute the maximum of additions in parallel. But on the other hand, if the batches are too small, the synchronization cost would increase, due to the bigger number of synchronization signals to manage. This is especially true as the granularity of doublings/halvings and additions (several hundreds of processor clock cycles) is too small in comparison with the cost of synchronization barriers and signals.

Q = k • P 1: Q ← O 2: for i from 0 to t -1 do 3: if ki = 1 then 4: Q ← Q + P 5:
The three following methods can be used to synchronize the two threads:

mutex. A mutex is a mutual exclusion lock provided by the pthread library used to synchronize threads. When a thread holds a mutex, another thread, trying to take it, is locked, waiting for the releasing of the mutex from the first thread. Mutexes are generally used to protect critical sections of code.

The cost of a lock or an unlock is about 150-200 processor clock cycles, which is almost negligible. signals: they are used in the inter-thread and inter-process communication.

A thread waiting for a signal is put in a sleeping state until another thread sends the corresponding signal. Then, the thread wakes up and goes on running. The sleeping state allows savings of resources which are then available for another process. In our experience and on our platform, the cost to send a signal is about 2000 clock cycles.

busy-waiting: this method consists in using a shared flag (in the global memory) and use it to keep the addition-thread in a busy-waiting loop while waiting for the producer-thread to output the next point and modify the flag. The main drawback of this method is to waste processor resources.

According to our experiments, signals are too costly compared to the two other techniques. The busy-waiting and mutex techniques almost give the same results in terms of performance, although the mutex method is slightly better in some cases. Thus we decided to use exclusively mutexes.

Proposed synchronization method. Our strategy was to avoid the use of mutex synchronization as much as possible. We chose to use only one single mutex: at the very beginning of the computation the mutex keeps the additionthread in an inactive state while a first batch of doublings or halvings is computed by the producer-thread. At the end of the computation of this batch, the producer-thread releases the mutex and pursues the whole sequence of doubling without performing any further locking on the mutex. This approach is depicted in Figure 1.

The correctness of the final result depends on the size of the first batch of points before the mutex releasing, which ensures that the writings of the point stored in shared memory by the doubling thread precedes the reading of the same point by the addition thread. If this batch is too small and in case of long sequence of zeros in the binary or NAF scalar representation, one can meet a violation of the read-after-write dependency, and the computation is not correct. To avoid this configuration, we carefully tuned this batch size in order to have the error rate as close as possible to zero. In our test results shown below, this error rate is limited to less than 1%. This is a compromise chosen in order to limit the first batch of doublings/halvings size, and to get the best performances. But at this step, such an error rate remains unacceptable. In order to eliminate these errors, we added a test on the addition-thread. In the producer-thread, we used a variable which is stored in global memory as the loop counter. This allows to check if the addition processed uses a point which has already been computed by the producer-thread, i.e. the read-afterwrite dependency is ensured. The cost of this test is almost negligible, although the use of a global memory counter is not totally free. When an error is detected (that is to say a read-after-write dependency violation), we break the additionthread loop, and launch a sequential computation of k • P . Due to the small error rate, the cost of this rescue computation, which frequency is near zero, is negligible on average. Algorithm 4 presents an algorithmic formulation in the case of Right-to-left Double-and-add scalar multiplication of this approach, including the elimination of the error computations due to a synchronization failure.

Impact of scalar recoding. In the sequential case, it is a useful technique to recode the scalar using NAF and W-NAF to speed-up the computation (as previously mentioned in Subsection 2.1 page 4). In the parallel algorithms, the situation is different. Indeed, the NAF and W-NAF recodings reduce the number of additions performed by the addition-thread. This fact can be seen when analyzing the amount of computations performed by the two threads. We can evaluate this amount using the results given Table 2 and Table 1 in the case of curves over E(F p ), and using the results given Table 2 and Table 3 in the case of curves over E(F 2 m ). For simplicity we assumed that S = 0.8M in F p and that a squaring and square root are negligible in F 2 m and that the cost of a quadratic solver is roughly 1M . The resulting complexities are given in Table 4. Table 4, we remark that, generally, the amount of computation of the additionthread is larger than the producer-thread for the binary coding. When using the NAF recoding the amount of computation of the two threads are roughly the same. Finally, the use of W-NAF makes the amount of computation of the addition-thread significantly smaller than the producer thread. This means that when using W-NAF recoding, the addition-thread progresses faster and even would have to wait for the producer-thread to output new points. But in any case, the addition-thread terminates after the producer-thread. Moreover in the W-NAF case, the post-computations delay the output of the results after the end of the producing process, since in the parallel algorithms, this final reconstruction cannot be done before the end of the parallelized additions.

These remarks are confirmed by the chronogram given in Figure 2 which shows the different timings required by each thread related to the recoding used for the execution of the parallelized halve-and-add for scalar multiplication in E(F 2 233 ). This fact leads us to opt for the NAF recoding for our implementations.

In addition to this choice, and in order to improve the performances, we implement a variable initial batch size of doublings/halvings in the producerthread. Indeed, the number of additions performed by the addition-thread depends on the Hamming weight of the NAF representation of the scalar (i.e. the non-zero digits). As stated previously, a read-after-write dependency violation can appear if this batch is too small. But if the Hamming weight of the scalar is higher, the risk of this dependency violation is lower, and the batch size can be reduced in this case. This improvement applies only when the addition-thread has roughly the same running time as the producer thread. This concerns the Double-and-add approach over F p and F 2 m but not the Halve-and-add approach (cf. Table 4).

When possible, we use the variable initial batch size in the producer-thread.

Four-thread parallel version over binary elliptic curve

Over binary field, the parallelization proposed by Taverne et al. in [START_REF] Taverne | Speeding Scalar Multiplication over Binary Elliptic Curves using the New Carry-Less Multiplication Instruction[END_REF] splits the scalar multiplication into two independent threads. Specifically, they split the t-bit scalar k = k 1 + k 2 where k 1 and k 2 are as follows

k = (k t 2 t-+ . . . + k ) k1 + (k -1 2 -1 + . . . + k 0 2 -) k2 .
(5) In general is close to t/2 and represents the length of the Halve-and-add subkey.

Then the computations can be parallelized into one thread computing k 1 P with the Double-and-add algorithm and a second thread computing k 2 P with the Halve-and-add algorithm.

We propose to combine the approach of Taverne et al. with the parallelization approach discussed in Subsection 3.2. This results in a four-thread algorithm: the partial scalar multiplication k 1 P is computed with the parallel two-thread algorithm // Double-and-add and k 1 P is computed with the parallel two-thread algorithm // Halve-and-add. This four-thread approach is shown in Figure 3. This approach increases the level of parallelization, but it also requires additional thread launching and management. Therefore, this algorithm works better on large fields, as it will be shown in the next section.

Recode k k is split in two subkeys (> 0 powers of 2, and ≤ 0 powers of 2). Our implementations of the four-thread of Figure 3 use the following strategies for thread launching and synchronization:

-The threads are launched in this order: 1) the Halve-and-add producerthread, which launches 2) the Double-and-add producer-thread which launches 3) the Halve-and-add addition-thread which finally launches 4) the Doubleand-add addition-thread. -The recoding of the scalar is done by 3) the Halve-and-add addition-thread before launching 4) the Double-and-add addition-thread. -Due to the delay of thread launching and key recoding computation, it is not necessary to use mutexes with initial batch size of halvings or doublings of points for each producer-thread. -In order to eliminate computation errors due to synchronization failure, we use the same method as the one described Subsection 3.2 for. Thus, we have two global memory counters, one for the doubling producer-thread and one for the halving producer-thread. Each addition-thread compare its own counter with the global counter of its corresponding, and launches a partial rescue computation if a synchronization dependency violation is detected.

Timings

The platform used for the experimentations is an Optiplex 990 DELL , with a Linux 12.04 operating system. The processor is an Intel Core i7 -2600 Sandy Bridge 3.4GHz, which owns four physical cores. The code is written in C language, compiled with gcc version 4.6.3. The Hyperthreading BIOS and also the Turbo-boost options have been deactivated on our platform in order to measure the performances as accurately as possible.

Since the operating system has the possibility to preempt the resources in order to launch another task, we avoid such difficulties by choosing to run our codes in a recovery mode shell. But we noticed that the codes generally run well in normal operating system conditions too, although perturbations may be observed in a few cases.

The performances are measured using one hundred batches of 2000 computations, each batch with a different random scalar. The minimum value of each batch is considered and the average value gives the performance. With this measurement process, we take into account the variations due to the different Hamming weights of the scalars.

For each case we provide the detailed duration of each thread. We notice that, generally, the overall computation finishes around several tens of thousands cycles after the producer-thread. These timings might correspond to the delayed start of the addition-thread (due to the initial batch size of points computed by the producer thread) and the synchronization and thread management time. For the four-thread versions, the given value corresponds to the scalar bit size of the Halve-and-add computation (cf. equation ( 5)). We have evaluated the overhead due to error management due to wrong error synchronization and it represents, in average, roughly 2-6% of the overall computation time.

Concerning the results over F 2 233 , we remark that the four-thread version is not competitive. This might be due to the synchronization and thread creation and management cost. Furthermore, the speed-ups with the two thread versions are not very important.

Over F 2 409 , the situation is different since the four-thread version is now better: it requires 324395 clock cycles whereas the two-thread parallel W-NAF Double/halve-and-add necessitates 347000 clock cycles (6.6% improvement). The speed-ups provided by the two-thread versions is also more important: between 15% (Double-and-add case) and 19.5% (Halve-and-add case).

Concerning the results over F p , we first notice that a scalar multiplication over a Jacobi Quartic is faster than over a Weierstrass curve. This corroborates the complexities of the curve operations shown in Table 1. We also notice that the tested two-thread parallelization provides performance improvements of around 15% to 17% compared to the NAF sequential Double-and-add approach.

Comparison. We give in Table 6 some published results in the literature. Over, F p , the work of Longa is on Intel Core 2 with p = 2 256 -189 and Hamburg is over a Sandy Bridge with p = 2 252 -2 232 -1 and has smaller key. The other works deal with the same processor and on the same fields as the one considered in this paper. We can see that, in the case of E(F 2 233 ) our two-thread approach is not competitive with best know results. In the cases of E(F 2 409 ), the proposed approach improves by 9.4% the previous best known timings reported in [START_REF] Taverne | Speeding Scalar Multiplication over Binary Elliptic Curves using the New Carry-Less Multiplication Instruction[END_REF]. Finally, the timing provided in [START_REF] Hamburg | Fast and compact elliptic-curve cryptography[END_REF] is better than the timing obtain by our method, but Hamburg uses a slightly smaller field and key size. On the other hand, we improve the best known results for curve defined for a 128 bit security level.

Conclusion

In this work, we have considered parallelized software implementations of scalar multiplication kP over E(F 2 m ) and E(F p ). We first have considered the parallelization suggested by Moreno et Hasan in [START_REF] Moreno | SPA-Resistant Binary Exponentiation with Optimal Execution Time[END_REF] which splits the right-to-left scalar multiplication into two threads: one producer-thread computing 2 i P or

8 :

 8 return (Q) Synchronization between threads. Both parallelization (right-to-left Doubleand-add, Algorithm 3 and Halve-and-add, Algorithm 2) are classical producerconsumer configurations.

Fig. 1 .

 1 Fig. 1. Synchronization and thread processing for our ECSM implementation

Fig. 2 .Algorithm 4

 24 Fig. 2. Chronogram of the halve-and-add computation with binary, NAF and W-NAF scalar representation over B233

Fig. 3 .

 3 Fig. 3. Four-thread algorithm.

Table 4 .

 4 Complexity of the two threads for a t-bit scalar coded in binary, NAF and W-NAF, in multiplication number.

		Double-and-add		Double-and-add	Halve-and-add
			over Fp			over F 2 m			over F 2 m
	Recoding	producer-addition-post-producer-addition-post-producer-addition-post-
		thread	thread comp.	thread	thread comp.	thread	thread comp.
	binary	6.2tM	5.1tM	0	4tM	6.5tM	0	2tM	4tM	0
	NAF	6.2tM	3.4tM	0	4tM	4.33tM	0	2tM	2.66tM	0
	W-NAF (w = 4)	6.2tM	2.04tM 33M	4tM	2.6tM	39M	2tM	1.6tM	39M

Acknowledgement: We would like to thank Christophe Nègre for his valuable and helpful comments.

This work has been suported by a PHD grant from PAVOIS project (ANR 12 BS02 002 01).

We give the parameters of the curves used in the experimentation in Appendix:

-Appendix A.1 for the B233 and B409 binary field curves; -Appendix A.2 for the Weierstrass prime field curve; -Appendix A.3 for the Jacobi-Quartic prime field curve.

Table 5 shows the results of the proposed parallel strategies for scalar multiplication implementations. The above timings include the error detection and correction due to erroneous thread synchronization.

Scalar

Curve Security processor Method Cycles multiplication

Hamburg [START_REF] Hamburg | Fast and compact elliptic-curve cryptography[END_REF] Montgomery 126 Intel core i7 SB Montgomery ladder 153000 Langley [START_REF] Langley | C25519 code[END_REF] Curve25519 128 Intel core i7 SB Montgomery ladder 229000 1 compiled and run on our platform.

Table 6. Performance comparison with the state of the art 2 -i P for i = 1, . . . , t and one addition-thread which accumulates these points to compute kP . We have proposed a lightweight approach for thread synchronization. In addition, in order to avoid remaining computation error due to dependency violation, we proposed a low cost checking method of the synchronization between threads with a rescue computation. We have also evaluated the best approach for the scalar recoding in this context. In the special case of E(F 2 m ) we have combined this approach to the parallelized Double/halve-andadd approach of [START_REF] Taverne | Speeding Scalar Multiplication over Binary Elliptic Curves using the New Carry-Less Multiplication Instruction[END_REF]. The experimental results show that these parallelization techniques provide some speed-up on elliptic curve scalar multiplication computations compared to previously best known implementations. Indeed, over prime field and binary fields, in most cases the parallelization provides an improvement of roughly 15% on the computation time.

A Appendix: Curve Parameters

A.1 Elliptic curves over binary field

The curve equation is:

Jean-Marc Robert

The parameters are for B233:

b = 0x00000066 647ede6c 332c7f8c 0923bb58 213b333b 20e9ce42 81fe115f 7d8f90ad, r = 0x00000100 00000000 00000000 00000000 0013e974 e72f8a69 22031d26 03cfe0d7.

where the order of the curve is n × h. For B409 we have:

b = 0x0021a5c2 c8ee9feb 5c4b9a75 3b7b476b 7fd6422e f1f3dd67 4761fa99 d6ac27c8 a9a197b2 72822f6c d57a55aa 4f50ae31 7b13545f, r = 0x01000000 00000000 00000000 00000000 00000000 00000000 000001e2 aad6a612 f33307be 5fa47c3c 9e052f83 8164cd37 d9a21173.

A.2 Weierstrass curve over prime field

The curve equation is:

The parameters are:

r is the prime order of P .

A.3 Jacobi Quartic curve over prime field

The curve equation is:

The parameters are:

The Jacobi Quartic curve is isomorphic to the following Weierstrass elliptic curve:

where: a = (-16 -3θ 2 )/4 and b = -θ 3 -aθ. Hence, in our case:


Elisa Rebolini
email: rebolini@lct.jussieu.fr

Julien Toulouse
email: julien.toulouse@upmc.fr

Andrew M Teale

Trygve Helgaker

Andreas Savin
email: savin@lct.jussieu.fr

Excitation energies along a range-separated adiabatic connection

We present a study of the variation of total energies and excitation energies along a rangeseparated adiabatic connection. This connection links the non-interacting Kohn-Sham electronic system to the physical interacting system by progressively switching on the electron-electron interactions whilst simultaneously adjusting a one-electron effective potential so as to keep the ground-state density constant. The interactions are introduced in a range-dependent manner, first introducing predominantly long-range, and then all-range, interactions as the physical system is approached, as opposed to the conventional adiabatic connection where the interactions are introduced by globally scaling the standard Coulomb interaction. Reference data are reported for the He and Be atoms and the H2 molecule, obtained by calculating the short-range effective potential at the full configuration-interaction level using Lieb's Legendre-transform approach. As the strength of the electron-electron interactions increases, the excitation energies, calculated for the partially interacting systems along the adiabatic connection, offer increasingly accurate approximations to the exact excitation energies. Importantly, the excitation energies calculated at an intermediate point of the adiabatic connection are much better approximations to the exact excitation energies than are the corresponding Kohn-Sham excitation energies. This is particularly evident in situations involving strong static correlation effects and states with multiple excitation character, such as the dissociating H2 molecule. These results highlight the utility of long-range interacting reference systems as a starting point for the calculation of excitation energies and are of interest for developing and analyzing practical approximate range-separated density-functional methodologies.

I. INTRODUCTION

Range-separated density-functional theory (see, e.g., Ref. 1) constitutes an interesting alternative to standard Kohn-Sham (KS) density-functional theory (DFT) [2,3]. In the standard KS approach, the physical interacting electronic Hamiltonian is replaced by an effective non-interacting Hamiltonian. By contrast, in range-separated DFT, the physical Hamiltonian is instead replaced by a partially interacting Hamiltonian that incorporates the long-range part of the electron-electron interaction. This corresponds to an intermediate point along a range-separated adiabatic connection [1,[START_REF] Savin | Recent Dev. Mod. Density Funct. Theory[END_REF][START_REF] Yang | [END_REF][6][7]. The KS Hamiltonian is linked to the physical Hamiltonian by progressively switching on the long-range part of the two-electron interaction, whilst simultaneously modifying the oneelectron potential so as to maintain a constant ground-state density. The ground-state energy of the physical system can then be extracted from the ground state of the long-range interacting Hamiltonian by using a short-range density functional describing the complementary short-range part of the electron-electron interaction. Be aware that this range-separated manner of introducing the interaction is not the usual way of performing the adiabatic connection, where the Coulomb interaction is instead scaled by a multiplicative constant going from 0 to 1.

Several short-range density-functional approximations have been developed [1,[START_REF] Savin | Recent Dev. Mod. Density Funct. Theory[END_REF][8][9][10][11][12][13] and a diverse range of approaches for calculating the ground state of the long-range interacting Hamiltonian have been explored. To aid in the description of static (or strong) correlation effects, which are poorly treated by standard density functionals, configuration-interaction [1,[START_REF] Savin | Recent Dev. Mod. Density Funct. Theory[END_REF]7,[14][START_REF] Savin | Recent Adv. Density Funct. Theory[END_REF][START_REF] Leininger | [END_REF][17], multiconfiguration self-consistent-field (MCSCF) [18][19][20], density-matrix functional theory (DMFT) [21][22][23], and constrained-pairing mean-field theory [24,25] descriptions of the long-range interacting systems have been employed. To treat van der Waals interactions, second-order perturbation theory [26][27][28][29][30][31][32][33][34][35][36][37], coupled-cluster theory [11,13,[38][39][40], and random-phase approximations [41][42][43][44][45][46][47][48][49][50][51] have been used successfully.

Electronic excitation energies can also be calcu-lated in range-separated DFT by using the linearresponse approach with a time-dependent generalization of the static ground-state theory [52]. In this case, the excitation energies of the long-range interacting Hamiltonian act as starting approximations that are then corrected using a shortrange density-functional kernel, just as the KS excitation energies act as starting approximations in linear-response time-dependent density-functional theory (TDDFT). Several such range-separated linear-response schemes have been developed, in which the short-range part is described by an approximate adiabatic semi-local density-functional kernel and the long-range linear-response part is treated at the Hartree-Fock [52][53][54][55], MCSCF [52,55], second-order polarization-propagator approximation (SOPPA) [55], or DMFT [56] levels. These schemes aim at overcoming the limitations of standard linear-response TDDFT applied with usual adiabatic semi-local approximations for describing systems with static correlation [57], double or multiple excitations [58], and Rydberg or chargetransfer excitations [59,60].

For the purpose of analyzing linear-response range-separated DFT approaches, it is desirable to have accurate reference values of the excitation energies of the long-range interacting Hamiltonian along the range-separated adiabatic connection [cf. Eq. ( 5)]. In this work, we provide and analyze reference data for the He and Be atoms and the H 2 molecule. The short-range oneelectron potentials required to keep the grounddensity constant along a range-separated adiabatic connection [cf. Eq. (6)] are determined at the full configuration-interaction (FCI) level using Lieb's Legendre-transform approach [61][62][63]. The excited-state energies of the long-range interacting Hamiltonian along the adiabatic connection [cf. Eq. (10)] are then calculated using the FCI method. Several accurate ground-state calculations have been performed in the past along the standard adiabatic connection [62][63][64][65][66][67] and rangeseparated adiabatic connections [1,6,[67][68][69] for small atomic and molecular systems, but accurate calculations of excited-state energies along adiabatic connections are very scarce-see, however, Refs. 62, 70.

The remainder of this paper is organized as follows. In Section II, range-separated DFT is briefly reviewed and the definition of the excited states along the range-separated adiabatic connection is introduced. In Section III, the behaviour of the excited-state energies near the two endpoints of the adiabatic connection, the KS system and the physical system, is studied analytically. After giving computational details in Section IV, results along the full adiabatic-connection path are presented and discussed in Section V. Finally, some concluding remarks are made in Section VI.

II. RANGE-SEPARATED DENSITY-FUNCTIONAL THEORY

In range-separated DFT (see, e.g., Ref. 1), the exact ground-state energy of an N -electron system is obtained by the following minimization over normalized multi-determinantal wave functions Ψ:

E 0 = min Ψ Ψ| T + Vne + Ŵ lr,µ ee |Ψ + Ēsr,µ Hxc [n Ψ ] . (1)
This expression contains the kinetic-energy operator T , the nuclear-electron interaction operator Vne = v ne (r)n(r)dr expressed in terms of the density operator n(r), and a long-range (lr) electron-electron interaction operator

Ŵ lr,µ ee = 1 2 w lr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 , (2)
expressed in terms of the pair-density operator n2 (r 1 , r 2 ). In the present work, we use the errorfunction interaction

w lr,µ ee (r 12 ) = erf(µr 12 ) r 12 , (3)
where µ controls the range of the separation, with 1/µ acting as a smooth cut-off radius. The corresponding complementary short-range (sr) Hartree-exchange-correlation density functional Ēsr,µ Hxc [n Ψ ] is evaluated at the density of Ψ: n Ψ (r) = Ψ|n(r)|Ψ .

The Euler-Lagrange equation for the minimization of Eq. (1) leads to the (self-consistent) eigenvalue equation

Ĥlr,µ |Ψ µ 0 = E µ 0 |Ψ µ 0 , (4)
where Ψ µ 0 and E µ 0 are the ground-state wave function and associated energy of the partially interacting Hamiltonian (with an explicit long-range electron-electron interaction) Ĥlr,µ = T + Vne + Ŵ lr,µ ee + V sr,µ Hxc .

(

) 5
It contains the short-range Hartree-exchangecorrelation potential operator, evaluated at the density n 0 (r) = Ψ µ 0 |n(r)|Ψ µ 0 , which is equal to the ground-state density of the physical system for all µ,

V sr,µ Hxc = vsr,µ Hxc [n 0 ](r)n(r)dr, (6)
where vsr,µ

Hxc [n](r) = δ Ēsr,µ Hxc [n] δn(r) . (7)
For µ = 0, Ĥlr,µ reduces to the standard noninteracting KS Hamiltonian, ĤKS , while for µ → ∞ it reduces to the physical Hamiltonian Ĥ: ĤKS = Ĥlr,µ=0 = T + Vne + VHxc , (8) Ĥ = Ĥlr,µ=∞ = T + Vne + Ŵee .

Varying the parameter µ between these two limits, Ĥlr,µ defines a range-separated adiabatic connection, linking the non-interacting KS system to the physical system with the ground-state density kept constant (provided that the exact short-range Hartree-exchange-correlation potential vsr,µ Hxc (r) is used).

In this work we also consider the excited-state wave functions and energies of the long-range interacting Hamiltonian

Ĥlr,µ |Ψ µ k = E µ k |Ψ µ k , (10)
where Ĥlr,µ is Hamiltonian in Eq. ( 5), with the short-range Hartree-exchange-correlation potential evaluated at the ground-state density n 0 . In range-separated DFT, these excited-state wave functions and energies provide a natural first approximation to the excited-state wave functions and energies of the physical system. For µ = 0, they reduce to the single-determinant eigenstates and associated energies of the non-interacting KS Hamiltonian,

ĤKS |Φ KS k = E KS k |Φ KS k , (11)
while, for µ → ∞, they reduce to the excitedstate wave functions and energies of the physical Hamiltonian

Ĥ|Ψ k = E k |Ψ k . (12)
We emphasize that, even with the exact (shortrange) Hartree-exchange-correlation potential, the total energies E KS k (E µ k ) are not the exact energies of the physical system but the total energies of a non-interacting (partially interacting) fictitious system of electrons with Hamiltonian ĤKS ( Ĥlr,µ ). Note also that, since the ionization energy is related to the asymptotic decay of the ground-state density, the ionization energy of the Hamiltonian in Eq. ( 10) is independent of µ and is equal to the ionization energy of the physical system. This is an appealing feature since it sets the correct energy window for bound excited states. Finally, note that the excitation energies ∆E µ k = E µ k -E µ 0 calculated from Eq. ( 10) constitute a starting point for range-separated linear-response theory based on the time-dependent generalization of Eq. ( 1) [52].

In this section, we study analytically the behavior of the excited-state energies E µ k as a function of the range-separation parameter µ close to the endpoints of the adiabatic connection: the KS system at µ = 0 and the physical system at µ → ∞. This study will aid in the understanding of the numerical results presented in Section V.

A. Excited-state energies near the Kohn-Sham system

We first derive the expansion of the excited-state energies near µ = 0, to see how the KS energies are affected by the introduction of the long-range electron-electron interaction. We assume that the system is spatially finite.

We rewrite the long-range interacting Hamiltonian of Eq. ( 5) as

Ĥlr,µ = ĤKS + Ŵ lr,µ ee -V lr,µ Hxc , (13)
with the long-range Hartree-exchange-correlation potential operator

V lr,µ Hxc = VHxc -V sr,µ Hxc = v lr,µ Hxc (r)n(r)dr. ( 14
)
The expansion of the long-range two-electron interaction is straightforward [1] (valid for µr 12 ≪ 1)

w lr,µ ee (r 12 ) = erf(µr 12 ) r 12 = 2µ √ π + µ 3 w lr, (3) ee
(r 12 ) + O(µ 5 ), (15)
with

w lr,(3) ee (r 12 ) = - 2 3 √ π r 2 12 . (16)
Note that the first term in the expansion of w lr,µ ee (r 12 ) in Eq. ( 15) is a spatial constant, 2µ/ √ π, which shows that what we call the long-range interaction does in fact contain also a contribution at short range [1]. Next, the expansion of the longrange Hartree-exchange-correlation potential

v lr,µ Hxc (r) = δE lr,µ Hxc [n] δn(r) (17)
can be determined from the expansion of the corresponding energy functional E lr,µ Hxc [n]. As derived in Ref. 1, the expansion of the Hartree-exchange part begins at first order and may be written as

E lr,µ Hx [n] = µ √ π n KS 2 (r 1 , r 2 )dr 1 dr 2 + µ 3 2 n KS 2 (r 1 , r 2 )w lr,(3) ee (r 12 )dr 1 dr 2 +O(µ 5 ). ( 18
)
where n KS 2 (r 1 , r 2 ) is the KS pair density, while the expansion of the correlation part only begins at sixth order (assuming a non-degenerate KS ground state)

E lr,µ c [n] = 0 + O(µ 6 ). ( 19
)
If the functional derivative of E lr,µ Hx [n] is taken with respect to density variations that preserve the number of electrons, δn(r)dr = 0, then the first-order term in Eq. ( 18) does not contribute due to the fixed normalization of the KS pair density, n KS 2 (r 1 , r 2 )dr 1 dr 2 = N (N -1). The derivative is then defined up to an additive (µ-dependent) constant C µ , which can be fixed by requiring that a distant electron experiences zero potential interaction in Eq. ( 13), amounting to setting the zeroenergy reference. The linear term in µ in the longrange Hartree-exchange-correlation potential can then be determined as follows.

To first order in µ, the long-range electronelectron interaction tends to a constant, 2µ/ √ π. A distant electron (with 1 ≪ r 12 ≪ 1/µ) then experiences a constant interaction 2(N -1)µ/ √ π with the remaining N -1 other electrons. This constant must be exactly compensated by the longrange Hartree-exchange-correlation potential in Eq. ( 13), so that its first-order term in µ must also be 2(N -1)µ/ √ π. The expansion of v lr,µ Hxc (r) therefore takes the form

v lr,µ Hxc (r) = 2(N -1)µ √ π + µ 3 v lr, (3)
Hxc (r) + O(µ 5 ), (20)
where v lr,(

Hxc (r) is the third-order contribution. Combining Eqs. [START_REF] Savin | Recent Adv. Density Funct. Theory[END_REF] and (20), we arrive at the following expansion of the long-range interacting Hamiltonian of Eq. ( 13):

Ĥlr,µ = ĤKS + µ Ĥlr,(1) + µ 3 Ĥlr,(3) + O(µ 5 ), (21) with a constant first-order correction Ĥlr,(1) = -

N (N -1) √ π (22)
and the following third-order correction

Ĥlr,(3) = Ŵ lr,(3) ee - V lr,(3) Hxc , (23)
Ŵ lr,(3

) ee = 1 2 w lr,(3) ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 , (24)
V lr,(3) Hxc = v lr,( 3)

Hxc (r)n(r)dr. ( 25
)
Since the first-order correction in the Hamiltonian is a constant, it does not affect the associated wave functions. The expansion of the wave functions therefore begins at third order in µ:

Ψ µ k = Φ KS k + µ 3 Ψ (3) k + O(µ 5 ). ( 26
)
Using normalized KS wave functions Φ KS k |Φ KS k = 1, the expansion of the total energy for the state k is then

E µ k = E KS k - N (N -1) √ π µ +µ 3 Φ KS k | Ĥlr,(3) |Φ KS k + O(µ 5 ). ( 27
)
The first-order contribution is the same for all states, cancelling out in the differences between the energies of two states. As a result, the corrections to the KS excitation energies are third order in µ.

For closed shells, the expansion of the difference between the singlet and triplet energies associated with the single excitation i → a can be obtained by applying Eq. ( 27) with the spin-adapted KS wave functions 1 Φ KS = Φ KS i→a + Φ KS ī→ā / √ 2, for the singlet state, and 3,1 Φ KS = Φ KS ī→a , for the triplet state with spin projection M S = 1. Only the two-electron term then contributes:

∆E µ,1-3 i→a = 2µ 3 ia| ŵlr,(3) ee |ai + O(µ 5 ) = 8µ 3 3 √ π | i|r|a | 2 + O(µ 5 ), (28)
where we have used r 2 12 = r 2 1 +r 2 2 -2r 1 •r 2 . The appearance of the transition dipole moment integral in Eq. (28) means that, for an atomic system, the singlet-triplet energy splitting appears at third order in µ if the difference between the angular moment of the orbitals ϕ i and ϕ a is ∆ℓ = +1 or -1. Otherwise, the splitting appears at a higher order in µ.

B. Excited-state energies near the physical system

We now derive the asymptotic expansion of the excited-state energies when µ → ∞, which shows how the exact excited-state energies are affected by the removal of the very short-range part of the electron-electron interaction.

For this purpose, we rewrite the long-range interacting Hamiltonian of Eq. ( 5) as

Ĥlr,µ = Ĥ -Ŵ sr,µ ee + V sr,µ Hxc , ( 29
)
where Ĥ is the Hamiltonian of the physical system, Ŵ sr,µ ee = 1 2 w sr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 (30) is the short-range electron-electron interaction operator defined with the complementary errorfunction interaction

w sr,µ ee (r 12 ) = erfc(µr 12 ) r 12 , (31)
and V sr,µ Hxc is the short-range Hartree-exchangecorrelation potential operator in Eq. ( 6). The first term in the asymptotic expansion of w sr,µ ee (r 12 ) can be written in terms of a delta function [1] (valid for µr 12 ≫ 1)

w sr,µ ee (r 12 ) = π µ 2 δ(r 12 ) + O 1 µ 3 , (32)
while the expansion of vsr,µ Hxc (r) = δ Ēsr,µ Hxc [n]/δn(r) can be obtained from that of Ēsr,µ

Hxc [n]. As derived in Ref. 1, the expansion of the long-range Hartreeexchange energy is

E sr,µ Hx [n] = π 2µ 2 n KS 2 (r, r)dr + O 1 µ 4 , ( 33
)
where n KS 2 (r, r) is the KS on-top pair density, while the expansion of the long-range correlation energy is

Ēsr,µ c [n] = π 2µ 2 n 2,c (r, r)dr + O 1 µ 3 , (34)
where n 2,c (r, r) is the on-top correlation pair density of the physical system. Therefore, the expansion of the short-range Hartree-exchangecorrelation potential takes the form vsr,µ

Hxc (r) = 1 µ 2 vsr,( -2)
Hxc (r) + O 1 µ 3 , (35)
where vsr,(-2) Hxc (r) is the µ -2 contribution formally obtained by taking the functional derivative of Eqs. (33) and (34).

Substituting Eqs. (32) and (35) into Eq. ( 29), we obtain the asymptotic expansion of the long-range interacting Hamiltonian as

Ĥlr,µ = Ĥ + 1 µ 2 Ĥlr,(-2) + O 1 µ 3 , (36)
where Ĥlr,(-2) = -Ŵ sr,(-2) ee

+ V sr,( -2)
Hxc is composed of an on-top two-electron term and a one-electron term:

Ŵ sr,(-2) ee = π 2 n2 (r, r)dr, (37)
V sr,(-2) Hxc = vsr,(-2)

Hxc (r)n(r)dr. ( 38
)
The expansion of the Hamiltonian in Eq. ( 36) suggests a similar expansion for the excited-state wave

functions, Ψ µ k = Ψ k + µ -2 Ψ (-2) k + O(µ -3
). However, as shown in Ref. 71, this expansion is not valid for r 12 ≪ 1/µ. The contribution of the wave function for small r 12 to the integral for the total energy E µ k = Ψ µ k | Ĥlr,µ |Ψ µ k nevertheless vanishes in the limit µ → ∞, and the asymptotic expansion of the total energy of the state k is

E µ k = E k + 1 µ 2 Ψ k | Ĥlr,(-2) |Ψ k + O 1 µ 3 , ( 39
)
where the wave function Ψ k is normalized to unity.

IV. COMPUTATIONAL DETAILS

Calculations have been performed for the He and Be atoms and for the H 2 molecule with a development version of the DALTON program [START_REF] Dalton | a molecular electronic structure program[END_REF][START_REF] Aidas | [END_REF], using the implementation described in Refs. 63 and 69. First, a FCI calculation was performed to determine the exact ground-state density within the basis set considered, followed by a Lieb optimization [62] of the short-range potential v sr,µ (r) = v ne (r) + vsr,µ Hxc (r) also at the FCI level to reproduce the FCI ground-state density in the presence of the long-range electron-electron interaction w lr,µ ee (r 12 ). The FCI excited-state energies were then calculated using the partially interacting Hamiltonian with the interaction w lr,µ ee (r 12 ) and effective potential v sr,µ (r).

The Lieb maximization was performed using the short-range analogue of the algorithm of Wu and Yang [74], in which the potential is expanded as

v sr,µ (r) = v ne (r) + v sr,µ ref (r) + t b t g t (r). ( 40
)
where the reference potential is the short-range analogue of the Fermi-Amaldi potential

v sr,µ ref (r) = N -1 N n 0 (r ′ )w sr ee (|r -r ′ |)dr ′ , (41)
calculated for a fixed N -electron density n 0 , to ensure the correct asymptotic behaviour. The same Gaussian basis set {g t } is used for the expansion of the potential and the molecular orbitals.

The coefficients b t are optimized by the Newton method, using a regularized Hessian with a truncated singular-value-decomposition cutoff of 10 -7 for He and 10 -6 for Be and H 2 .

Even-tempered Kaufmann basis sets [75] and uncontracted correlation consistent Dunning basis sets [76] augmented with diffuse functions were tested extensively for the He atom, especially to converge the lowest P state. No significant differences were observed using the two basis sets and Calculations were performed for about 30 values of µ between 0 to 10 bohr -1 (with about half the points between 0 and 1 where the energies vary the most). Cubic spline interpolation has been used on this calculated data when plotting the total and excitation energies as a function of µ. For later use, analytical expressions were also fitted to the calculated total energies and excitation energies. The forms used in the fitting were chosen to satisfy the expansions at small and large µ values as presented in Eqs. (27) and (39). The details of these fits are given in the supplementary material [77].

V. RESULTS AND DISCUSSION

A. Helium atom

The total energies of the ground state 1 1 S and of the first Rydberg-like singlet and triplet S and P excited states of the He atom are plotted as a function of the range-separation parameter µ in Figure 1. At µ = 0, the KS non-interacting total energies are obtained. Being sums of orbital energies with a resulting double counting of electron repulsion, these quantities are well above the µ in bohr -1 Excitation energies in hartree total energies of the physical system (higher by about 1 hartree). When the long-range electronelectron interaction is added by increasing µ from µ = 0, the total energies decrease linearly with µ with a slope of -2/ √ π, in accordance with the linear term in the expansion of Eq. ( 27) for N = 2. For larger µ values, the total energy curves flatten and approach the energies of the physical system asymptotically as 1/µ 2 as µ → ∞, in accordance with Eq. (39). The total energies along the adiabatic connection are poor approximations to the total energies of the physical system unless the range-separation parameter µ is large. Specifically, µ 6 is required to be within 10 mhartree of the exact total energies.

1 1 S → 2 3 S 1 1 S → 2 1 S 1 1 S → 1 3 P 1 1 S → 1 1 P Figure 2: Excitation energies ∆E µ k = E µ k -E µ 0 (
The lowest singlet and triplet excitation energies are plotted in Figure 2. The KS singlet and triplet excitation energies are degenerate and, as already observed for a few atomic systems in Refs. 78-80, are bracketed by the singlet and triplet excitation energies of the physical system. As µ increases from µ = 0, the excitation energies vary as µ 3 since the linear term in Eq. ( 27) cancels out for energy differences. The singlet-triplet degeneracy is lifted and the excitation energies converge to the exact singlet and triplet excitation energies when µ → ∞. Whereas a monotonic variation of the excitation energy with µ can be observed for the singlet and triplet 1S→ 2S excitations and for the triplet 1 1 S → 1 3 P excitation, a non-monotonic variation is observed for the singlet 1 1 S → 1 1 P excitation. This behaviour could be an artefact of the basis-set expansions (either orbital or potential), noting that a similar behaviour was ob- served for other excitations in a smaller basis set and was removed by enlarging the basis set (the basis set dependence of the singlet 1 1 S → 1 1 P excitation energy is given in the supplementary material [77]). In line with previous observations in Refs. [START_REF] Umrigar | Electronic Density Functional Theory[END_REF]80] for the KS system, the excitation energies for Rydberg-type states along the adiabatic connection are rather good approximations to the excitation energies of the physical system (the maximal error is about 0.02 hartree at µ = 0 for the triplet 1 1 S → 2 3 S excitation), becoming better and better for high-lying states as they must eventually converge to the exact ionization energy.

The singlet-triplet energy splittings for the 2S and 1P states are plotted in Figure 3. The expansion at small µ of Eq. ( 28) predicts the singlettriplet splitting to increase as µ 3 for the 1P state since it corresponds to the 1s → 2p excitation in the KS system, so that ∆ℓ = 1. By contrast, the singlet-triplet splitting should increase at most as µ 5 for the 2S state since it corresponds to the 1s → 2s excitation in the KS system, so that ∆ℓ = 0. This difference is clearly visible in Figure 3, where the 2S curve for the singlet-triplet splitting initially increases more slowly than the 1P curve.

B. Beryllium atom

The total energies of the ground state 1 1 S and of the valence singlet and triplet 1P excited states of the Be atom are plotted in Figure 4. The KS total energies are approximately 6 hartree above the physical energies. At small µ, an initial slope of -12/ √ π is observed for all states, in accordance with Eq. ( 27) with N = 4. However, convergence to the physical energies with increasing µ is much slower than for the He atom, owing to the short inter-electronic distances in the Be 1s core region, which are consequently probed at larger µ values.

The singlet and triplet excitation energies are plotted in Figure 5. As for He, the KS excitation energies are bracketed by the singlet and triplet excitation energies of the physical system. Not surprisingly, the KS excitation energies are poorer approximations to the exact excitation energies for these valence excitations in Be than for the Rydberg excitations in He. As µ increases, the KS excitation energies rapidly converge to the physical excitation energies. Clearly, the slow convergence of the core energies does not affect the convergence of the valence excitation energies.

Close to the KS system, at µ = 0, the excitation energies are quite sensitive to the introduction of a small portion of electron-electron interaction in the Hamiltonian, which may be interpreted as a sign of static correlation. For µ ≈ 0.4 -0.5, a typical µ value in range-separated DFT calculations [18,81], the calculated excitation energies are significantly better approximations to the exact excitation energies than are the KS excitation energies. This observation justifies range-separated multi-determinantal linear-response DFT calculations, which take these excitation energies as a starting point.

C. Hydrogen molecule

The first few excitation energies of H 2 at the equilibrium bond distance are plotted against µ in Figure 6. As for the atoms, the valence excitations energies vary much more along the adiabatic connection than do the Rydberg-like excitation energies. Note also that the energetic ordering of the states changes along the adiabatic connection. With our choice of basis set, we also observe that the higher singlet excitation energies do not depend monotonically on µ, approaching the physical limits from above, as observed for He. Again, the excitation energies around µ ≈ 0.4 -0.5 represent better approximations to the exact excitation energies than the KS excitation energies.

Finally, we consider the interesting case of the dissociation of the H 2 molecule. The first excitation energies at three times the equilibrium distance are shown in Figure 7. With increasing bond distance, the 1σ g and 1σ u molecular orbitals become degenerate. Consequently, the KS excitation energy for the single excitation 1σ g → 1σ u goes to zero. Moreover, the KS excitation energy for the double excitation (1σ g ) 2 → (1σ u ) 2 also goes to zero (albeit more slowly). This behaviour is in contrast to that of the physical system, where only the excitation energy to the triplet 1 3 Σ + u state goes to zero, whilst those to the singlet 1 1 Σ + u state and the 2 1 Σ + g state (the latter connected to the double excitation in the KS system) go to finite values.

Clearly, the excitation energies of KS theory are poor approximations to the exact excitation energies, making it difficult to recover from these poor starting values in practical linear-response µ in bohr -1 Excitation energies in hartree µ in bohr -1 Excitation energies in hartree TDDFT calculations. As µ increases from µ = 0, the excitation energies to the singlet 1 1 Σ + u and 2 1 Σ + g states vary abruptly, rapidly approaching the physical values. This sensitivity to the inclusion of the electron-electron interaction is a clear signature of strong static correlation effects, emphasizing the importance of a multi-determinantal description in such situations. At µ ≈ 0.4-0.5, the 1 1 Σ + u and 2 1 Σ + g excitation energies, although still too low, are much better approximations than the KS excitation energies, constituting a strong motivation for range-separated multi-determinantal approaches in linear-response theory.

1 1 Σ + g → 1 3 Σ + u 1 1 Σ + g → 1 1 Σ + u 1 1 Σ + g → 2 3 Σ + g 1 1 Σ + g → 2 1 Σ + g 1 1 Σ + g → 1 3 Πu 1 1 Σ + g → 1 1 Πu Figure 6: Excitation energies ∆E µ k = E µ k -E µ 0 (
1 1 Σ + g → 1 3 Σ + u 1 1 Σ + g → 1 1 Σ + u 1 1 Σ + g → 2 1 Σ + g

VI. CONCLUSIONS

We have studied the variation of total energies and excitation energies along a range-separated adiabatic connection, linking the non-interacting KS system (µ = 0) to the physical system (µ → ∞) by progressively switching on the long-range part of the electron-electron interaction with the range-separation parameter µ, whilst keeping the ground-state density constant. This behaviour is of interest for the development and analysis of range-separated DFT schemes for the calculation of excitation energies, such as the linear-response range-separated schemes of Refs. 52, 53, 55. Reference calculations were performed for the He and Be atoms and the H 2 molecule. Except when µ is large, the ground-and excited-state total energies along the adiabatic connection are poor approximations to the corresponding energies of the physical system. On the other hand, the excitation energies are good approximations to the excitation energies of the physical system for most of the adiabatic connection curve, except close to the KS system (µ = 0). In particular, the excitation energies obtained at µ ≈ 0.4 -0.5, typically used in range-separated DFT calculations, are significantly better approximations to the exact excitation energies than are the KS excitation energies. This behaviour appears to be particularly evident for situations involving strong static correlation effects and double excitations, as observed for the dissociating H 2 molecule.

These observations suggest that the excitation energies of the long-range interacting Hamiltonian in range-separated DFT may be useful as first estimates of the excitation energies of the physical system. However, if one cannot afford to use large µ values (µ > 2 -3), these excitation energies should be considered only as starting approximations, suitable for correction by, for example, linear-response range-separated theory.

In future work, we will utilize the present reference data to assess the approximations made in practical linear-response range-separated schemes, where the long-range contribution is treated, for example, at the Hartree-Fock, MCSCF or SOPPA levels of theory, while the short-range part is described by semi-local density-functional approximations. We will also use the results of this work to guide the development of timeindependent range-separated DFT methods for the calculation of excitation energies as alternatives to linear-response schemes-in particular, for methods based on perturbation theories [START_REF] Filippi | [END_REF]82] or extrapolations [83,84] along the adiabatic connection.

Figure 1 :

 1 Figure 1: Ground-and excited-state total energies E µ k (in hartree) of the He atom as a function of µ (in bohr -1 ). The total energies of the physical system E k = E µ→∞ k are plotted as horizontal dotted lines. The slope at µ = 0 is shown by the black dashed line for the first excited state.

  Figure 2: Excitation energies ∆E µ k = E µ k -E µ 0 (in hartree) of the He atom as a function of µ (in bohr -1 ). The excitation energies of the physical system ∆E k = ∆E µ→∞ k

Figure 3 :

 3 Figure 3: Singlet-triplet energy splittings (in hartree) for the He atom as a function of µ (in bohr -1 ).

Figure 4 :µ 1 S → 1 1 PFigure 5 :

 4115 Figure 4: Ground-and excited-state total energies E µ k (in hartree) of the Be atom as a function of µ (in bohr -1 ). The total energies of the physical system E k = E µ→∞ k are plotted as horizontal dotted lines. The slope at µ = 0 is shown in dashed line.

  Figure 6: Excitation energies ∆E µ k = E µ k -E µ 0 (in hartree) of the H2 molecule at the equilibrium internuclear distance as a function of µ (in bohr -1 ). The excitation energies of the physical system ∆E k = ∆E µ→∞ k

Figure 7 :

 7 Figure 7: Excitation energies ∆E µ k = E µ k -E µ 0 (in hartree) of the H2 molecule at 3 times the equilibrium internuclear distance as a function of µ (in bohr -1 ). The excitation energies of the physical system ∆E k = ∆E µ→∞ k

Acknowledgments E. R. and J. T. gratefully acknowledge the hospitality of the Centre for Theoretical and Computational Chemistry (CTCC), University of Oslo, where part of this research was done. E. R. also thanks A. Borgoo and S. Kvaal for helpful discussions. T. H. acknowledges support from the Norwegian Research Council through the CoE Centre for Theoretical and Computational Chemistry (CTCC) Grant No. 179568/V30 and the Grant No. 171185/V30 and through the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683. A. M. T. is grateful for support from the Royal Society University Research Fellowship scheme.