Oscar Nierstrasz

Marcus Denker

Lukas Renggli

Model-Centric, Context-Aware Software Adaptation ⋆

Software must be constantly adapted to changing requirements. The time scale, abstraction level and granularity of adaptations may vary from short-term, fine-grained adaptation to long-term, coarsegrained evolution. Fine-grained, dynamic and context-dependent adaptations can be particularly difficult to realize in long-lived, large-scale software systems. We argue that, in order to effectively and efficiently deploy such changes, adaptive applications must be built on an infrastructure that is not just model-driven, but is both model-centric and contextaware. Specifically, this means that high-level, causally-connected models of the application and the software infrastructure itself should be available at run-time, and that changes may need to be scoped to the run-time execution context. We first review the dimensions of software adaptation and evolution, and then we show how model-centric design can address the adaptation needs of a variety of applications that span these dimensions. We demonstrate through concrete examples how model-centric and context-aware designs work at the level of application interface, programming language and runtime. We then propose a research agenda for a model-centric development environment that supports dynamic software adaptation and evolution.

Introduction

It is well-known that real software systems must change to maintain their value [START_REF] Lehman | Program Evolution: Processes of Software Change[END_REF]. It is therefore curious to observe that the technology we use to develop software systems tends to hinder and inhibit change rather than to enable and support it [START_REF] Nierstrasz | Change-enabled software systems[END_REF]. Statically typed languages, for example, are based on the assumption that first-class values have fixed types that will not change, especially at run-time. Few mechanisms are available to developers to deal with the fact that interfaces do change over time, and real software systems may need to cope with different versions of the same libraries, possibly depending on the runtime context. Design patterns offer further evidence of ungainly workarounds that developers need to regain flexibility at run-time, for example to change the apparent behaviour of objects as a consequence of a change in state [START_REF] Gamma | Design Patterns: Elements of Reusable Object-Oriented Software[END_REF].

Long-lived, software intensive systems [START_REF]Report of the Beyond the Horizon thematic group 6 on Software Intensive Systems[END_REF] cannot always be modified in a static way. Furthermore, although certain kinds of anticipated adaptations can be built in by design as run-time configuration parameters, there are many kinds of dynamic adaptation that cannot be anticipated so easily. One canonical example of such an adaptation is run-time instrumentation: certain kinds of anomalies only manifest themselves with deployed software systems. As it is not possible to anticipate for all cases what and where to trace to observe the problematic behaviour, it may be necessary to dynamically adapt the running system. Other examples exist (such as adding new features to an always-running system), but the key characteristics remain the same -the software may need to be adapted dynamically, in a fine-grained way, while taking care not to disturb existing behaviour.

There are many important dimensions of software change. Let us just consider three of these that pose challenges for software development:

Timescale -Software is changed not only at the coarse scale of versions and releases, but also at a medium scale (e.g., start-up configuration) and at a fine scale (run-time adaptation and instrumentation). Particularly at the dynamic end, little support is available to developers aside from certain design patterns and relatively low-level reflective mechanisms. Granularity -Here too we see that software is changed not only at the coarse granularity of subsystems and packages, or the medium granularity of classes and methods, but also at a finer granularity within methods and procedures. Fined-grained, run-time adaptation of software must typically be anticipated by design, and necessitates the use of boilerplate code (e.g., case-based reasoning over anticipated scenarios) or design patterns (e.g., State or Strategy patterns). Unanticipated run-time adaptation will typically entail low-level techniques such as bytecode transformation. Scope -Changes may be globally visible, they may be localized to individual users, or they may depend on an even finer context. The same software entities may need to behave differently as the run-time context changes. Mobile applications, for example, may need to switch to a fall-back behaviour as services become unavailable. Run-time instrumentation of software entities, as another example, may need to be dynamically adapted if the same entities are used by the instrumentation layer itself (i.e., to avoid endless instrumentation loops) [START_REF] Denker | The meta in meta-object architectures[END_REF].

Although model-driven and round-trip engineering techniques have proved to be effective in maintaining the connection between high-level and low-level views of software systems, they do not especially address the problem of dynamic adaptation. We argue that it is necessary to go a step further from model-driven towards model-centric software, in which high-level, causally connected views of software and their application domain are available at run-time. In this paper we show several examples of run-time adaptation at the level of source code, so the "high-level models" appropriate to these applications take the form of ASTs that reflect the structure of software to be adapted. Furthermore, such systems must be context-aware in order to control the scope of adaptations and changes. In our examples we show how context can play an important role in software adaptation to control the scope of change. We argue that current programming technology offers only very weak support for developing context-aware applications, and that new research is urgently needed into novel context-oriented programming mechanisms [START_REF] Hirschfeld | Context-oriented programming[END_REF].

In this paper we make our case for model-centric, context-aware software adaptation by presenting two examples of platforms that adopt this approach. We show how the presence of sufficiently high-level models at run-time can enable very dynamic forms of context-dependent software adaptation.

In Section 2 we present Reflectivity, a relatively mature platform for dynamic, model-centric software adaptation. We have used Reflectivity extensively in various projects to support different forms of adaptation, such as run-time instrumentation, dynamic aspects, and software transactional memory. Next, in Section 3, we present ongoing work on Diesel, a lightweight language workbench which can be used to adapt the programming environment to support the expression of high-level application concepts by introducing numerous, lightweight domain-specific languages. We discuss further applications of these ideas and our vision for a research agenda in Section 4 and provide an overview of related work in Section 5. We conclude in Section 6 with some remarks on future work.

Reflectivity -a platform for model-centric software adaptation

In this section we present Reflectivity, a platform that supports dynamic adaptation of software by means of causally connected, high-level models of the source code [START_REF] Denker | Sub-method Structural and Behavioral Reflection[END_REF]. The purpose of this section is (i) to motivate the need for dynamic software adaptation for various applications such as runtime instrumentation, dynamic aspects, and software transactional memory, (ii) to motivate the need for better mechanisms to support context-dependent adaptation, and (iii) to demonstrate that sufficiently high-level models available at run-time (in this case ASTs causally connected to bytecode) facilitate run-time adaptation. Reflectivity is built on top of Smalltalk, since it already provides extensive support for run-time reflection, albeit at a relatively low-level of abstraction [START_REF] Denker | Sub-method Structural and Behavioral Reflection[END_REF]. Furthermore, Smalltalk provides full access to the implementation of its infrastructure, making it ideal for extensive experimentation. Any other language that supports run-time structural and behavioural reflection and access to the infrastructure would also be suitable.

A model for dynamic software adaptation

The particular challenge we are focusing on is support for dynamic, fine-grained and possibly context-dependent software adaption. Let us consider the canonical example of run-time instrumentation:

-We may need to install the instrumentation code dynamically in the running system because the phenomena we wish to study only occur in the deployed system (say, a web service). -The adaptation is fine-grained because we wish to monitor only part of a given method (say, conditional access to an authorization service). -The adaptation is context-dependent because we are only interested in monitoring calls made from a specific application, not others.

Other plausible scenarios, such as adding features to a running system, would serve as well for establishing our requirements.

In order to dynamically adapt software, we need a model to reason about it. In our run-time instrumentation scenario the following properties are important:

Abstraction Level: This model should be high-level, reflecting the language concepts we wish to instrument, rather than, say, the generated bytecode. Completeness: The model should represent the complete software, from coarsegrained structures like classes, methods down to sub-method structures such as variable accesses and method calls.

Although these properties may seem obvious, in most cases the representations used for software adaptation today do not satisfy them. The representation used is often plain (source) text. Modern development environments do better: here the code is represented with dedicated data-structures that better support code presentation (e.g. pretty printing) or code change (e.g. refactoring). But these data structures are those of the development environment, not of the language itself. They are not available to support run-time adaptation.

Runtime representations are often tailored solely towards execution, such as bytecode representations for Java or Smalltalk. Representations based on bytecode are low-level, and therefore suffer from a semantic mismatch with the core language concepts.

The reflective representation of the structure of software available in many modern object-oriented languages provides a high-level model for packages, classes and methods, but it lacks any representation of sub-method structure.

As we are especially interested in adaptation at runtime and by the system itself, we conclude that the model needs to have the following properties:

Self representation: The model of the software needs to be available from within the running system itself. Causal connection: When we change this model (either from the outside or from within the system), the behavior of the program needs to change. Conversely, when the system changes, the representation needs to change, too. The program needs to stay in sync with the model at all times. Meta-annotations: We need to be able to extend the representation to use it in many contexts and annotate it with meta-data. For example, different tools that deal with the structure of the system need slightly different information.

Annotations allow the programmer to associate meta-data with any node, making the existing AST-based representation extensible.

Models that have the properties of self representation and causal connection are called reflective. There is a long history of reflection in programming languages in general and in object-oriented languages in particular [START_REF] Tanter | Reflection and open implementations[END_REF]. Submethod reflection [START_REF] Denker | Submethod reflection[END_REF][START_REF] Denker | Sub-method Structural and Behavioral Reflection[END_REF] provides a model of the software that exhibits all the properties discussed above. Software is represented down to the statement level by a causally connected, annotated Abstract Syntax Tree (AST). Changes to the AST will be (lazily) propagated to the generated bytecode using runtime justin-time compilation. Semantics are given to annotations by an open compiler infrastructure: annotations are interpreted by dedicated compiler plug-ins.

We will now illustrate this approach by a series of examples.

Run-time instrumentation

With a causally connected model of the whole software, we can provide a framework to support instrumentation at runtime. The model chosen is that of partial behavioral reflection [START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF][START_REF] Denker | Higher abstractions for dynamic analysis[END_REF][START_REF] Röthlisberger | Unanticipated partial behavioral reflection: Adapting applications at runtime[END_REF]. The central notion is the Link. The link is set as an annotation to one or more nodes of our AST. The link points to a meta-object and can be parameterized to indicate which information is passed to the meta-object. In addition, we can set a condition to specify when a link is active and to specify if the meta-object is called before, after or instead of the original instruction. Figure 1 shows the interaction of the AST, the link and the meta-object. Links are specified as annotations on the AST. A compiler plugin transforms the AST before execution to take the links into account. A link thus results in code to be inserted in the program at the nodes where it is installed.

This model provides some interesting characteristics: it is completely dynamic due to just-in-time compilation at runtime. We can create links at runtime, configure them and install them in the system. Links can even be installed by other links, or they can remove or install themselves.

A side-effect of using the higher-level representation provided by the AST is improved performance. It is actually easier to generate more efficient code using the AST [START_REF] Denker | Submethod reflection[END_REF]. Furthermore, we only have to generate bytecode for those parts of the system that take part in the execution. The actual set of classes used by an application is generally smaller than the overall code base by an order of magnitude. As a result, dynamic code generation provides an additional performance benefit [START_REF] Denker | Sub-method Structural and Behavioral Reflection[END_REF].

Localization: annotating structure

We have seen both a structural model of our system, and a framework for dynamic behavioral reflection based on annotating the structural model with links. Now we will see how to manipulate this structure using behavioural reflection.

To make this possible, we need to be able to reference the structural model from the behavioral world. The simplest way to do this is to allow the nodes of the structural model to be meta-objects, as shown in Figure 2. This allows a node to be annotated before it is executed. This way one can easily realize tracing or feature analysis [START_REF] Denker | Supporting feature analysis with runtime annotations[END_REF]. For example, a simple code-coverage tool can be realized by installing a link on each AST node of interest. We just provide a method markExecuted to mark AST nodes as having been executed. The links are activated when the AST nodes are executed, and simply invoke this method to record the fact. link := GPLink new metaObject: #node; selector: #markExecuted.

Listing 1.1. Code-coverage analyzer realized with Reflectivity.

When we install the link on the node representing methods, we obtain method level coverage. But with our sub-method model, we can go a level deeper and even install the link on all assignments.

To improve performance, it is even possible to remove the tagging-link at runtime just after tagging the node. In this way, the method would, at the next execution, be recompiled to only call markExecuted on those nodes that have not yet been executed before. We can also take advantage of activation conditions, for example, to only tag nodes that are executed in the context of a unit-test.

The possibility of both installing and removing links at runtime allows for just-in-time annotation: links are installed on-demand on all methods that are to be executed next. This way an annotation can spread itself through the system, driven by the flow of execution itself. Examples like these make Reflectivity especially suitable for building self-monitoring and self-evolving systems.

Scoping the effect of changes

An interesting problem arises when instrumenting basic system classes. The instrumented code itself is used at runtime by the meta-object, leading to endless loops. This makes any use of reflection on basic system classes like Number or Array impractical. This problem can be seen with all reflective systems -a well-known example is CLOS [START_REF] Chiba | Avoiding confusion in metacircularity: The meta-helix[END_REF].

To solve this problem, we provide the possibility to scope the activation of links towards meta-level execution [START_REF] Denker | The meta in meta-object architectures[END_REF]. Links are parameterized with the level for which they are activated. This way we can restrict the introduced change towards, for example, base-level program execution. Note that the same mechanisms can be used to reason about execution of meta-level programs. For example, it may be the case that a profiler realized as a meta-object needs to be analyzed to improve performance. By restricting the link that activated the profiler to the meta-level, we can use the profiler on itself without the danger of endless loops.

Implementing higher-level dynamic language features

With partial behavioral reflection, we can easily adapt the programming language to support new language features. Very deep changes can be realized that normally require changes at the level of language implementation (i.e. the virtual machine).

Dynamic aspects. Partial behavioral reflection can serve as an efficient technique for implementing Aspect-Oriented Programming [START_REF] Rodríguez | Supporting dynamic crosscutting with partial behavioral reflection: a case study[END_REF]. We have used the dynamic features of Reflectivity to implement dynamic aspects, which are not woven into the code at compile time, but instead can be introduced and retracted at runtime. Compared with traditional runtime AOP implementation techniques, we can see some improvements. We can generate better code than typical bytecode transformation based approaches as we can leverage the higherlevel AST representation [START_REF] Denker | Sub-method Structural and Behavioral Reflection[END_REF]. We can leverage the link-conditions to efficiently control aspect activation at runtime.

Transactional memory. We have realized software transactional memory for a dynamic language [START_REF] Renggli | Transactional memory in a dynamic language[END_REF]. It is notable that this realization was done without any changes to the underlying virtual machine. With the help of Reflectivity we introduced a transactional execution context and were able to reify all state access to be handled by the transactional model implemented in the host language.

3 Diesel -an engine for bringing models closer to code

We have seen that it is important to have higher level models of software available at run-time in order to enable dynamic software adaptations. But what about the application logic itself? Models of the software structures can be far removed from the application domain. What we need to enable dynamic adaptation of application logic is to make application models more explicit in the code. One way of doing this is to raise the level of programming by means of specialized domain specific languages (DSLs).

This raises new, extrinsic problems, as a new DSL will not be able to automatically benefit from existing development tools available for the general-purpose host language. We need to be able to adapt the host language and environment to support new DSLs. As before, the adaptions must be fine grained because DSL code may be interspersed with regular source code, dynamic because we want to be able to change host compilers and tools on the fly, and context-dependent because DSL code may be restricted to certain parts of an application.

Diesel is a lightweight language workbench that closely integrates with the host language [START_REF] Fowler | Language workbenches: The killer-app for domainspecific languages[END_REF]. This enables developers to incrementally bend the syntax and semantics of the host language to suit their exact needs for a particular problem domain. As such, Diesel is an environment for developing domain specific languages (DSLs), with the aim of giving application developers more suitable abstractions than the host language provides. Contrary to other approaches, Diesel reuses the traditional compiler toolchain and closely integrates with the existing tools of the programming environment, such as editors, debuggers, inspectors, etc. A close integration is crucial to keep the abstraction gained through new language features. While language developers might want to toggle between a view on the original source and the transformed result, domain developers would like to stay at the abstraction level of their code at all times.

The host language of our implementation is Smalltalk, which provides us with a uniform development environment. The abstract code representation of the compiler is reused in all parts of the system and can thus take advantage of our extensions. For example, if we change the syntax in a specified part of the system, syntax highlighting and debugger continue to work. Application developers do not have to learn new tools, but continue to use the existing ones even if they mix multiple languages.

Example: modelling relationships

A common challenge in transforming UML models to code is how to implement relationships between objects. The problem has long been solved in relational databases [START_REF] Atkinson | Types and persistence in database programming languages[END_REF][START_REF] Nixon | Implementation of a compiler for a semantic data model: Experiences with taxis[END_REF], but none of today's mainstream languages provide first-class relationships [START_REF] Nelson | First class relationships for OO languages[END_REF]. If done by hand, it is easy to introduce subtle bugs that might be very hard to detect. With code generation, huge chunks of boilerplate code may appear that are hard to understand and impossible to change. In either case, debugging and maintaining the code is cumbersome, because developers not only have to think in terms of the high-level DSL code that they specify, but also in terms of the code that is generated.

We present an approach to solve this problem with Diesel. The example only handles 1:1 relationships, however it could easily be extended to support arbitrary n:m relationships.

To implement the write accessor next: of a double linked list, a developer or a code generator would write something like this:

1 Link>>next: aLink At run-time lines 2 and 4 are need to ensure that the inverse relations are properly updated. The actual assignment only occurs on line 3. Most parts of the code are not interesting to developers. It is an unnecessarily complex code fragment specifying an implementation detail of 1:1 relationships that is probably used at several places throughout the application.

We now replace the complex write accessor from above with a plain write accessor that does not update the opposite relationships:

Link>>next: aLink next := aLink
However we put high-level annotations next to the instance variable declaration of the Link class, to tell Diesel that all write access to these variables requires their respective inverse relationships to be updated: When compiling the method, Diesel will automatically generate the necessary boilerplate code around it. This generated code is never visible, not even in the debugger, to ensure that the developer can concentrate on the high-level model. The magic behind the transformation comes from a rule that has been added to the Diesel engine. Whenever source-code is parsed, translated and annotated these rules are processed to enable interaction with the compiler: The rule given above matches all parse tree nodes that assign an expression to an instance variable (line 2). The remaining code checks if the instance variable is annotated with an opposite annotation (lines 5-6) and then defines the transformation programmatically. This is done using a quasi-quoting mechanism [START_REF] Bawden | Quasiquotation in Lisp[END_REF] to build and inject the AST nodes that update the opposite relationship into the tree. Subsequently the tree is transformed to bytecodes, by the standard Smalltalk compiler. The handwritten and the transformed code result in identical bytecodes, therefore both approaches perform equally at runtime. In practice the minimal increase in compilation time due to the additional transformations can be neglected.

Scoping the effect of changes

In the above example the scope of the transformation is given at the level of parse tree nodes. Often such transformation rules only apply to a carefully chosen part of the system however. For example, the above transformation should be used in model code only, but not in the UI implementation. Diesel supports a wide variety of additional constraints that can be composed and added to transformation rules: packages, namespaces, classes, class hierarchies or even specific methods.

Furthermore arbitrary conditions can be added to the transformed parse-tree nodes, so that the transformation is only in effect if a certain runtime condition is met. The generated code can resort to the reflective capabilities of the system [START_REF] Rivard | Smalltalk: a reflective language[END_REF] and select the appropriate behaviour depending on the runtime context [START_REF] Costanza | Language constructs for context-oriented programming: An overview of ContextL[END_REF].

Towards a research agenda

Full reflection (i.e., with run-time intercession) has been widely available in dynamic programming languages for many years, notably in Smalltalk [START_REF] Rivard | Smalltalk: a reflective language[END_REF] and CLOS [START_REF] Kiczales | The Art of the Metaobject Protocol[END_REF]. Nevertheless, reflection has been commonly considered to be either too dangerous or too difficult to use for common programming tasks. Static languages such as Java and C++ offer a weaker form of reflection that only supports introspection at run-time (i.e., the possibility to examine but not to affect the model elements).

There is increasing pressure to adapt software systems at run-time. If the host programming language does not offer reflective features, programmers can be forced to adopt workarounds, such as reifying and interpreting model elements within their programs. This obviously places a heavy burden on developers who must build up this infrastructure themselves, possibly in ad hoc ways.

The Reflectivity framework simplifies the development of reflective applications by offering ASTs as relatively high-level, causally connected, run-time models of program elements. Reflection at the sub-method level is enabled since the deep structure of programs is captured, unlike in other approaches which stop at the method level. We have used the Reflectivity framework extensively for various purposes, including the analysis of software features [START_REF] Denker | Supporting feature analysis with runtime annotations[END_REF], dynamic monitoring of software entities from the IDE [START_REF] Röthlisberger | Exploiting runtime information in the IDE[END_REF], and even for the implementation of pluggable types, where type expressions are encoded as source code annotations and interpreted by compiler plug-ins [START_REF] Haldimann | Practical, pluggable types for a dynamic language[END_REF]. Despite these successes, Reflectivity is focused on reifying model elements of the host programming language, not those of the application domain. (One could also say that the application domain of Reflectivity is the host programming language.) What we are missing is Reflectivity for domain models. We envision a system where end users can change their domain models on the fly, without having to touch the host programming language [START_REF] Renggli | Magritte -a meta-driven approach to empower developers and end users[END_REF].

Traditionally domain concepts are translated and encoded in the source code in such a way that makes it difficult to reason about these concepts once the software is deployed. It is often difficult, for example, to find the software components responsible for a given end-user feature in the source code. In a model-driven approach, one would express domain concepts at the level of meta-models and models, and then generate code from these descriptions. In a straightforward approach, this still has the consequence that the domain models are no longer directly expressed in the code. If features must be dynamically adapted, there is no easy way to manipulate these model elements from the running system. Instead of generating code from models, we feel that it is necessary to bring models closer to code. This means that domain concepts should be expressed directly in the source code, rather than being encoded using concepts of the solution domain. In essence, rather than applications being model-driven, with models merely being used to generate code, we believe that they should be model-centric, with models being first-class entities that can be manipulated at run-time.

One way of bringing models closer to code is to provide higher-level, domain specific languages for model concepts more directly. One downside of this approach is the potential for the proliferation of DSLs, each with their own obscure syntax. DSLs should therefore be simple and lightweight. Another important downside is that the existing development environment will need to be adapted to work with each new DSL. Diesel addresses these problems by offering a lightweight framework for specifying simple DSLs that are transformed into the ASTs used by Reflectivity. The transformations are used to keep the development tools in sync with each DSL, so that editors and debuggers, for example, can present developers with the original DSL code rather than the generated host language code. In a sense, the model is the code.

Returning to the theme of software adaptation, we note that any adaptation manifests itself as a kind of software change, which is possibly intrusive and possibly context-dependent. As an example, consider software that should adapt its policies for ensuring changing requirements (e.g., related to concurrency control, or security, or transactional behaviour) dynamically according to (i) the run-time context (e.g., the presence of competing applications), and (ii) availability of related services (e.g., for optimistic or pessimistic transaction support). Such an application has clearly defined requirements but can only partially anticipate its run-time context and the nature of services available to meet those requirements.

Depending on the context, the same software will need to behave differently, and may need to be adapted in unanticipated ways.

We summarize the research directions that we see as essential to support dynamic software adaptation as follows:

Bring models closer to code. In order to enable dynamic software adaptation, models should be first-class, high-level artifacts available at run-time for both introspection and intercession. Structured source code and run-time annotations offer one light-weight technique to embed domain knowledge in source code [START_REF] Marschall | Persephone: Taking Smalltalk reflection to the sub-method level[END_REF]. Lightweight DSLs are another promising technique. Model-centric development. Rather than seeking ways to embed models in source code, perhaps we should replace the source code as an artifact and directly program with models. After all, third generation languages were once seen as a way to "generate (machine) code" from high-level specifications. Nowadays we consider programs in these languages to be the source code. We need to take the next step and jettison our third generation objectoriented languages in favour of models as source code. Environments to support model-centric development would concentrate on directly manipulating first class models and their meta-models. Models would be available at runtime to support the same kinds of adaptations available to the developer at development time. Context-oriented programming. COP refers to programming language mechanisms and techniques to support dynamic adaptation to context [START_REF] Hirschfeld | Context-oriented programming[END_REF]. Although many present-day applications need to be context-aware, contextdependent behaviour is generally programmed in ad hoc ways, due to the lack of support in modern programming languages. Some COP languages have been proposed [START_REF] Costanza | Language constructs for context-oriented programming: An overview of ContextL[END_REF][START_REF] González | Highly dynamic behaviour adaptability through prototypes with subjective multimethods[END_REF], and both Reflectivity and Diesel are examples of a frameworks that provide some degree of COP support, but research is very young, and there is no consensus how best to support COP in programming languages.

Related work

There is a long history of reflective programming languages, ranging from dynamic languages such as Lisp, Smalltalk, Scheme and CLOS, to static languages like C++ and Java, which provide a more limited form of run-time introspection rather than full intercession. All of these approaches are limited to models of the code base, and do not take models of requirements, design decisions or architecture into consideration. Over the years various approaches have attempted to keep high-level knowledge about software in sync with the software itself. The earliest examples of these is probably Literate Programming [START_REF] Knuth | Literate Programming. Stanford, California: Center for the Study of Language and Information[END_REF], in which documentation and source code are freely interspersed and maintained together.

In some cases Architectural Description Language (ADL) specifications are considered to be part of the running software system, rather than simply a higher-level description of it, as it is the case with Darwin [START_REF] Magee | Specifying distributed software architectures[END_REF]. Although ADLs provide a high-level interface for specifying and configuring components at an architectural level, there is not really any explicit representation of a model that is developed in tandem with the rest of the software.

Generative programming approaches [START_REF] Czarnecki | Generative programming: methods, tools, and applications[END_REF] produce software from higher-level descriptions using such mechanisms as generic classes, templates, aspects and components. A general shortcoming of these approaches is that the transformation is uni-directional -there is no way to go from the code back to higher-level descriptions. Round-trip engineering refers to approaches in which transformations are bi-directional [START_REF] Antkiewicz | Framework-specific modeling languages with round-trip engineering[END_REF]. Models and code are still considered to be separate artifacts, so models are not available at run-time for making adaptive decisions.

Case tools and 4GLs represent an attempt to simplify the generation and adaption of an application. However the main focus of these approaches was to generate code and not to consider models as executable artifacts of software development.

Model-driven engineering (MDE) [START_REF] Bézivin | Towards a precise definition of the OMG/MDA framework[END_REF][START_REF] Douglas | Guest editor's introduction: Model-driven engineering[END_REF] refers to a more recent trend in which application development is driven by the development of models at various levels of abstraction. Platform-independent models are transformed to platformspecific models, and eventually to code which runs on a specific platform. Generally these transformations are performed off-line, so models are not necessarily available to the run-time system, though some approaches support this [START_REF] Haustein | A model-driven runtime environment for web applications[END_REF].

The Eclipse Modeling Framework (EMF) [START_REF] Budinsky | Eclipse Modeling Framework[END_REF] provides facilities for manipulating models and generating Java source code from these models. Here too the focus is on models of source code, rather than on other views of a software system. The model and the code are still separate entities.

Naked objects [START_REF] Pawson | Naked Objects[END_REF] is an approach to software development in which domain objects and software entities are unified. Business logic is encapsulated in the domain objects and the user interface is completely generated from these domain objects. In this approach the domain model and the executing runtime are tightly coupled. Although naked objects address the earlier complaint against approaches which separate domain models from the source code or the running system, they do not offer any help in integrating other views of the software as it is being developed (i.e. requirements models, architectural views, and so on).

Aspect-Oriented Programming (AOP) [START_REF] Kiczales | Aspect-Oriented Programming[END_REF] provides a general model for modularising cross cutting concerns. Join points define points in the execution of a program that trigger the execution of additional cross-cutting code called advice. Join points can be defined on the runtime model (i.e., dependent on control flow). Although AOP works at a sub-method level, it does not provide a structural model of the system or any other reflective capabilities. The goal of AOP is to modularize crosscutting concerns, not to provide a model for dynamic software adaptation.

Reflex [START_REF] Tanter | Partial behavioral reflection: Spatial and temporal selection of reification[END_REF] pioneered partial behavioral reflection in the context of Java. Here links are associated with so called hooksets, abstractions of operations at the bytecode level. Therefore, the structural model of Reflex is that of bytecode, not the higher-level AST. In addition, Reflex does not support meta-annotations on bytecode.

Context-Oriented Programming (COP) [START_REF] Hirschfeld | Context-oriented programming[END_REF] refers to programming language support for developing applications whose behaviour depends on the run-time context. Present prototypes of COP languages focus on mechanisms for adapting behaviour to context, but provide little support for reasoning about context at the model level.

Changeboxes [START_REF] Denker | Encapsulating and exploiting change with Changeboxes[END_REF] provide a mechanism to control the scope of change in a running system. Deployment and development versions of a running software system can co-exist without interfering. Mechanisms for merging differences and resolving conflicts, however, must be handled in an ad hoc fashion, as no fully general approach exists for all usage scenarios. Changeboxes currently operate at the level of source code changes. There is no notion of changes to higher-level models.

Unlike general-purpose programming languages, domain specific languages (DSLs) tend to be compact languages that provide appropriate notations and abstractions for a particular problem domain. It was shown that DSLs increase productivity and maintainability for specialized tasks [START_REF] Van Deursen | Little languages: Little maintenance?[END_REF]. DSLs are often categorized as being either homogenous (internal), where the host-language and the DSL are one and the same, or heterogeneous (external), where the two languages are distinct [START_REF] Sheard | Accomplishments and research challenges in meta-programming[END_REF]. Techniques have been proposed to define language and semantics for new DSLs [START_REF] Holger Krahn | Integrated definition of abstract and concrete syntax for textual languages[END_REF]. The idea of designing languages that embrace the addition of new DSLs has been a focus of research in the past [START_REF] Odersky | Scala language secification v. 2.4[END_REF][START_REF] Warth | OMeta: an object-oriented language for pattern matching[END_REF][START_REF] Tratt | Domain specific language implementation via compile-time metaprogramming[END_REF]. However integrating new languages into existing tools has been largely neglected.

In the 1990s there was considerable interest in the development of architectural description languages (ADLs) [START_REF] Shaw | Software Architecture: Perspectives on an Emerging Discipline[END_REF] to capture and express architectural knowledge of a software system. ADLs can be viewed as DSLs for describing the architecture of complex software systems. Many DSLs formalize architecture in terms of components, connectors, and the rules governing their composition [START_REF] Shaw | Software Architecture: Perspectives on an Emerging Discipline[END_REF]. This idea is also implicitly contained in the notion of scripting languages, which can be seen as DSLs for composing applications from components written in another, usually lower-level programming language [START_REF] Ousterhout | Scripting: Higher level programming for the 21st century[END_REF]. Despite this, the interplay between conventional object-oriented languages, ADLs, scripting languages and DSLs has not yet been thoroughly studied nor has it been exploited in practice.

Concluding remarks

Software systems are under increasing pressure to support run-time adaptation for localization, mobile platforms, dynamic service availability, and countless other context-dependent applications. Unfortunately mainstream programming languages and development environments focus on limiting and restricting change rather than enabling it. Specifically, dynamic, fine-grained and contextdependent software adaptation is not well-supported by modern development technology.

In this paper we have presented two ongoing research projects that illustrate the principle of model-centric, context-aware software adaptation, and we have outlined a number of promising research directions for further exploration.

Reflectivity is a mature, model-centric framework for dynamic software adaptation. Software structures are represented at run-time by means of abstract syntax trees, which enables both dynamic and fine-grained adaptation. Furthermore, by means of partial behavioural reflection, adaptations can be scoped to the dynamic context. As a typical example, run-time instrumentation can be scoped to the context of a given feature, so it is possible to identify which software components support a given feature.

Diesel is a new framework for bringing models closer to code by supporting the definition of lightweight DSLs that are transformed to the host programming language using the same high-level source code models as provided by Reflectivity. The development tools, such as editors and debuggers, are aware of the transformations, so developers can continue to work with the high-level models rather than the transformed code, even while debugging. Here too, runtime models are used to support fine-grained adaptations that can be scoped to specific parts of the application requiring these DSLs.

Reflectivity and Diesel both build on top of the same rich programming environment and therefore could profit from each other in the future. On one hand, the causal connection of model and code would help Diesel to automatically propagate language changes to all its users. Furthermore Diesel could apply language transformations that normally happen at compile-time from the dynamic world. On the other hand, Reflectivity could profit from a richer language infrastructure and the compiler plugins could reuse the Diesel transformations.

These are only two examples of promising research directions to support dynamic change and software adaptation. We argue that more research is needed to close the gap between models and code. In general, software systems need to be change-enabled -instead of limiting and restricting change, they should actively enable change by treating change as a first-class entity [START_REF] Nierstrasz | Analyzing, capturing and taming software change[END_REF][START_REF] Nierstrasz | Change-enabled software systems[END_REF]. Ultimately we want model-centric and context-oriented programming environments where we can directly manipulate models both during development time and run-time, and software can be dynamically adapted by context.

Fig. 1 .

 1 Fig. 1. The link-meta-object model

Fig. 2 .

 2 Fig. 2. Bridging structural and behavioral model

2

  next isNil ifFalse: [ next instVarNamed: 'prev' put: nil ].

3 next := aLink. 4 next

 34 isNil ifFalse: [ next instVarNamed: 'prev' put: self ]

  Link instanceVariables: #( next <opposite: prev> prev <opposite: next> )

  context at: ' variable'.

5 opposite:

 5 = variable annotationNamed: 'opposite'.

8

  ifFalse: [ ,variable instVarNamed: ,opposite put: nil ]).

9

  context addNodeAfter: ( ,variable isNil 10 ifFalse: [ ,variable instVarNamed: ,opposite put: self ]) ]

Acknowledgments. We gratefully acknowledge the financial support of the Swiss National Science Foundation for the project "Bringing Models Closer to Code" (SNF Project No. 200020-121594, Oct. 2008 -Sept. 2010).

We thank Tudor Gîrba and Jorge Ressia for their careful reviews of various drafts. We also thank the anonymous reviewers for their numerous suggestions on how to improve the presentation of this paper.


Claude-Pierre Jeannerod

Jingyan Jourdan-Lu

Christophe Monat

Guillaume Revy

How to square floats accurately and efficiently on the ST231 integer processor

Keywords: squaring, binary floating-point arithmetic, correct rounding, IEEE 754, instruction level parallelism, C software implementation, VLIW integer processor

We consider the problem of computing IEEE floatingpoint squares by means of integer arithmetic. We show how the specific properties of squaring can be exploited in order to design and implement algorithms that have much lower latency than those for general multiplication, while still guaranteeing correct rounding. Our algorithm descriptions are parameterized by the floating-point format, aim at high instruction-level parallelism (ILP) exposure, and cover all rounding modes. We show further that their C implementation for the binary32 format yields efficient codes for targets like the ST231 VLIW integer processor from STMicroelectronics, with a latency at least 1.75x smaller than that of general multiplication in the same context.

Introduction

The STMicroelectronics ST231 is a 32-bit, 64-register, 4-issue, embedded integer VLIW processor. It is currently used in consumer electronics devices (set-top-boxes, sound processors, camera systems...) mainly performing intensive media and signal processing. As this processor has integeronly register file and ALUs, all the single precision floatingpoint support is available through software emulation, based on the highly optimized FLIP 1.0 library [START_REF] Jeannerod | FLIP 1.0: a fast floating-point library for integer processors[END_REF]. This comprehensive support is efficient enough to allow application developers to use out-of-the box floating-point code, instead of converting it to fixed-point representations.

However, the implementation of floating-point operations on integer-only processors can lead to sub-optimal use of the computational resources: for example, all binary op-erators f (x, y) called with the same argument x = y lead to redundant unpacking of the floating-point format, and to useless checks for special cases, and fail to exploit some properties of f (x, x). Since the compiler can detect all static cases of such occurrences, this leads to the idea of specializing operators to gain additional performance.

In this paper we focus on the specialization of the multiplication operator f : (x, y) → x × y into a square operator x → x 2 . Squares of floating-point values are ubiquitous in scientific computing and signal processing, since they are intensively used in any algorithm requiring the computation of Euclidean norms, powers, sample variances, etc. [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF][START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF].

We give a thorough study of how to design efficient software for IEEE floating-point squaring on targets like the ST231, that is, by means of integer arithmetic and logic only, and with high ILP exposure.

Our first contribution is to show how the specific properties of squaring can be exploited in order to refine the IEEE specification of multiplication, to deduce definitions of special input and generic input that are suitable for implementation, and to optimize the generic path and the special path for latency. This analysis is done for all rounding modes and presented in a parameterized fashion, in terms of the precision and the exponent range of the input/output floating-point format.

Second, this analysis allows us to produce a complete portable C code for the binary32 format and each rounding mode. On the ST231 processor, the result is a latency of 12 cycles for rounding 'to nearest even,' which is 1.75x faster than the 21 cycle latency of the multiply operator of FLIP 1.0; for the other rounding modes, the speedups are even higher and range from 1.9 to 2.3. Also, the average number of instructions issued every cycle lies between 3.4 and 3.5, thus indicating heavy use of the 4 issues available.

Third, we report on some experiments involving non-IEEE variants and square-intensive applications. We show that relaxing the IEEE requirements (finite math only, no support of subnormals) saves at most 1 cycle in the context of the ST231. We also show that for applications like the Euclidean norm, the sample variance, and binary powering the practical impact of our fast squarer reflects well the best that can theoretically be achieved.

For squaring in integer / fixed-point arithmetic several optimized hardware designs have been proposed; see for example [3, §4.9] as well as [START_REF] Walters | Combined unsigned and two's complement hybrid squarers[END_REF][START_REF] Gök | Integer squarers with overflow detection[END_REF] and the references therein. However, for squaring in IEEE floating-point arithmetic much less seems to be available and to the best of our knowledge, no optimized design has been presented and analyzed in the details as we do here, be it in hardware or in software. Furthermore, the implementation of squaring for the binary32 format and the ST231 processor outlined in [START_REF] Raina | FLIP: a Floating-point Library for Integer Processors[END_REF] does not support subnormal numbers, is available only for rounding 'to nearest even,' and has a latency of 27 cycles, which is 2.25x more than our 12 cycle latency.

Notation. For a real number r we write ⌊r⌋ for the largest integer not greater than r, and ⌈r⌉ for the smallest integer not less than r. If r satisfies further 0 ≤ r < 2 then its binary expansion i≥0 r i 2 -i with r i ∈ {0, 1} for all i will be simply written (r 0 .r 1 r 2 . . .) 2 . For a floating-point datum x, we write |x| for the absolute value of x, assuming implicitly that x is not NaN. Finally, we shall often use the following square bracket notation: for any true-or-false statement S, let [S] be 1 if S is true, 0 otherwise.

Outline. This paper is organized as follows. After some reminders on the ST231 processor and IEEE binary floating point in Section 2, we show in Section 3 how to specialize to squaring the IEEE specification of multiplication, deduce suitable definitions of generic and special input, and give a high-level algorithmic view of the squaring operator. Sections 4 and 5 then detail our algorithms for handling, respectively, generic and special input by means of integer arithmetic, and give the corresponding C implementation for the binary32 format. The performances of this implementation on the ST231 processor are reported in Section 6, and we conclude in Section 7. The detailed proofs of Properties 1 to 9 have been postponed to Appendix A for more legibility.

Background

STMicroelectronics' ST231 processor

The targeted processor is the ST231 from STMicroelectronics. Since it is a scoreboarded VLIW, there is no dynamic instruction dispatch contrary to what could achieve an equivalent 4-way superscalar architecture: reaching the best performance relies heavily on compiler efficiency.

Instruction scheduling is done after code selection, that must be aggressive to favor high ILP. Specifically, the if-conversion optimization that replaces conditional branches by conditional 'slct' instructions helps creating large straight-line code regions that are subject to efficient scheduling. The latency of this 'slct' instruction is 1 cycle.

Up to four independent instructions can be bundled into a 'syllable,' achieving the maximum throughput of four instructions per cycle. One constraint is that the immediate value that can be encoded as part of an instruction is limited to 9-bit signed: a larger 32-bit constant, named extended immediate, can be built by consuming an additional instruction syllable in the bundle. This reduces the actual parallelism available locally, but this compares favorably with other mechanisms such as loading the constant from memory, that incurs a performance impact on the data cache side of the machine.

Only a very small subset of the ST231 instruction set is needed by our squaring algorithms: logical and bitwise operators (&&, ||, &, |), relational operators (<=, >=, >), bitwise shift operators (<<, >>), addition, subtraction, maximum, and minimum of (un)signed integers (+, -, max, maxu, minu), and finally a multiply operator mul giving the higher half ⌊AB/2 32 ⌋ of the exact product of two 32-bit unsigned integers A and B. Except for mul whose latency is 3 cycles, the latency of each of these operators is 1 cycle.

IEEE 754-2008 binary floating point

Binary floating-point data. As specified in the IEEE 754-2008 standard [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], binary floating-point data consist of quiet or signaling not-a-numbers (qNaN, sNaN), signed zeros (+0, -0), signed infinities (+∞, -∞), and finite nonzero binary floating-point numbers having the following form: given a precision p and an exponent range [e min , e max ],

x = (-1) s • m • 2 e , ( 1a
)
where s is either 0 or 1, and where m = (m 0 .m 1 . . . m p-1 ) 2 and e min ≤ e ≤ e max . (1b)

Any such number must in fact be either normal (m 0 = 1) or subnormal (m 0 = 0 and e = e min ). Thus, writing α for the smallest positive subnormal number and Ω for the largest normal number, we have

α = 2 emin-p+1 and Ω = (2 -2 1-p ) • 2 emax .
Assumptions on e max , e min , and p. We shall assume that e max = 2 w-1 -1 for some positive integer w, (2a)

e min = 1 -e max , 2 ≤ p < e max , (2b)
and we shall write k = w + p.

All the binaryk formats of the IEEE 754-2008 standard satisfy these assumptions; see [7, Table 3.5] and, for a proof that 2 ≤ p < e max , see Appendix A. This allows us to carry out all the analysis in a parameterized way and specialize later to a given format, for example the binary32 format: w = 8, e min = -126, e max = 127, p = 24.

Standard encoding into k-bit integers. For the binaryk format the standard encoding of x in (1) is via a k-bit unsigned integer X whose bit string satisfies

X = [s|E w-1 • • • E 0 |m 1 • • • m p-1 ], E = w-1 i=0 E i 2 i , (3)
with E = e -e min + m 0 . Zeros, infinities, and NaNs are encoded via special values of X. In particular, writing |X| = X mod 2 k-1 , we have

x =          +0 iff X = 0, +∞ iff X = 2 k-1 -2 p-1 , sNaN iff 2 k-1 -2 p-1 < |X| < 2 k-1 -2 p-2 , qNaN iff |X| ≥ 2 k-1 -2 p-2 . ( 4
)
Correct rounding. Besides floating-point data and their encoding into integers, the standard [7, §4.3] defines rounding modes • to map any real number y to a unique floatingpoint datum x = •(y). In radix 2 four rounding modes are required: to nearest even (RN), down (RD), up (RU), and to zero (RZ). In the case of the square operator we shall restrict to RN, RD, and RU, since RZ(y) = RD(y) when y ≥ 0.

Exceptions and flags. Finally, the standard defines five exceptional situations (invalid operation, division by zero, overflow, underflow, inexact) and requires that they shall be signaled by raising some status flags. For square, all exceptions but 'division by zero' can occur, just like for multiply. However, since the status flags are currently not set in the multiply operator of the FLIP library [START_REF] Jeannerod | FLIP 1.0: a fast floating-point library for integer processors[END_REF] to which we compare, we have not implemented them for square either.

Specification and high-level algorithm

Squaring being a special case of general multiplication x × y, it is fully specified by the IEEE 754-2008 standard: given a rounding mode • and assuming x = y, the result r prescribed by [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] for x × y is as follows:

r =      |x| if x ∈ {±0, ±∞}, qNaN if x is NaN, •(x 2 ) otherwise. ( 5
)
This specification shows that r is essentially known in advance when x is zero, infinity, or NaN. However, in the particular case of squaring, if |x| is nonzero but "small enough" then the rounded value •(x 2 ) will always be equal to a tiny constant. Similarly, if |x| is finite but "large enough" then •(x 2 ) will always be equal to a huge constant. The rest of this section studies such properties of •(x 2 ) in order to refine the specification (5) and deduce a high-level algorithm.

Let min • and max • denote, respectively, the minimum value and maximum value of •(x 2 ) for |x| in [α, Ω]. Using the monotonicity, for x ≥ 0, of the map x → x 2 together with the definitions of rounding and overflow in [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], one may check that these values are as follows:

• RN RD RU min • +0 +0 α max • +∞ Ω +∞ (6)
For which values of x are such extremal values attained? To answer this, let us define the following two quantities

α ′ = 2 ⌊(emin-p)/2⌋ and Ω ′ = 2 (emax+1)/2 . ( 7
)
We give below three properties regarding these quantities.

Property 1. The values α ′ and Ω ′ defined in [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] are normal floating-point numbers such that α ′ < Ω ′ .

Property 2. For • ∈ {RN, RD, RU} and x a finite nonzero floating-point number, one has

•(x 2 ) = max • iff |x| ≥ Ω ′ .
The interval [Ω ′ , Ω] thus defines the widest input range on which max • is achieved. Interestingly, this range is the same for all our rounding modes, which makes things simpler from the implementer point of view. Note also that Ω ′ is an integer power of two because of (2a). For min • the situation is slightly more complex, as this minimum value is achieved on an input range [α, α ′ • ] whose upper bound now depends on •. However, the property below shows that one can suppress this dependency by restricting to input ranges whose upper bound has the form 2 i for some integer i. Property 3. The value α ′ in ( 7) is the largest integer power of two such that, for every finite nonzero floating-point number

x in [α, α ′ ), •(x 2 ) = min • for • ∈ {RN, RD, RU}.
The main outcome of Properties 2 and 3 is the following specification of floating-point squaring, which refines (5):

r =                    +0 if x = ±0, min • if α ≤ |x| < α ′ , •(x 2 ) if α ′ ≤ |x| < Ω ′ , max • if Ω ′ ≤ |x| ≤ Ω, +∞ if x = ±∞, qNaN if x is NaN. (8)
This brings a natural distinction between two kinds of input:

Definition 1. Input x is called generic if α ′ ≤ |x| < Ω ′ ,

and special otherwise.

By Property 1 every subnormal input is special, so that generic input consist of normal numbers only. The corresponding output •(x 2 ) must be finite because of the "only if" part in Property 2, but it can be (sub)normal or zero.

Let us now define

C spec = [x is special]. ( 9
)
This condition allows us to translate (8) into the following high-level algorithmic description, which shows that squaring essentially reduces to three independent tasks T i :

evaluate the condition C spec in (9) [T 1 ] if (C spec ) then handle special input as in (8) [T 2 ] else compute •(x 2 ) [T 3 ]
Thanks to if-conversion, the generated assembly for the above algorithm will consist of a straight-line piece of code computing the result R i of each task T i and ending with a 'slct' instruction that selects

R 2 if R 1 is true, R 3 otherwise.
For the design and implementation of each task we shall proceed in two steps as in [START_REF] Bertin | Techniques and tools for implementing IEEE 754 floating-point arithmetic on VLIW integer processors[END_REF]: assuming unbounded parallelism we optimize the a priori most expensive task first, namely task T 3 (see Section 4), and then only T 1 and T 2 , by trying to reuse as much as possible what was computed for T 3 (see Section 5). The latency of 'slct' being of 1 cycle, the lowest latency we can expect for squaring thus is 1 cycle more than that of T 3 .

Computing correctly-rounded squares

In this section we consider the computation of •(x 2 ) for x generic, that is, x as in [START_REF] Bertin | Techniques and tools for implementing IEEE 754 floating-point arithmetic on VLIW integer processors[END_REF] and such that

α ′ ≤ |x| < Ω ′ . (10)
By Property 1 such an x is normal, and thus 1 ≤ m < 2.

Parameterized formula for •(x 2 )

Normalized representation of the exact square. A first step towards the computation of

•(x 2 ) consists in normal- izing the representation m 2 • 2 2e of x 2 implied by (1a). Let c = m ≥ √ 2 . ( 11
)
Defining

m ′ = m 2 • 2 -c and e ′ = c + 2e then yields the unique pair (m ′ , e ′ ) ∈ R × Z such that 1 ≤ m ′ < 2 and x 2 = m ′ • 2 e ′ . ( 12
)
This is the so-called normalized representation of the exact square. Tight bounds for e ′ are given by the next result. These bounds for e ′ are the best possible ones:

• The lower bound is attained when |x| = α ′ . It is equal to e min -p-ǫ with ǫ = p is odd . One has ǫ = 0 for the standard binary32 format, and ǫ = 1 for the standard binary16, binary64, and binary128 formats [7, §3.6]. • The upper bound e max is attained for example when x = (1.1) 2 •2 (emax-1)/2 , which satisfies [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] and whose square is (1.001) 2 • 2 emax . Also, the tight lower bound on e ′ is less than e min for p ≥ 2, so that both situations e ′ ≥ e min and e ′ < e min do occur.

Correctly-rounded value of the exact square. When e ′ ≥ e min the normalized representation [START_REF] Raina | FLIP: a Floating-point Library for Integer Processors[END_REF] already allows to express the correctly-rounded value •(x 2 ). In this case x 2 lies in the range [2 emin , 2 emax+1 ) and, since m ′ ∈ [1, 2),

•(x 2 ) = •(m ′ ) • 2 e ′ = •(m 2 • 2 -c ) • 2 c+2e . ( 13a
)
When e ′ < e min the exact square ranges in (0, 2 emin ). In this case, we first set the exponent to e min and only then round the resulting scaled significand in fixed point:

•(x 2 ) = •(m ′ • 2 -(emin-e ′ ) ) • 2 emin = •(m 2 • 2 -(emin-2e) ) • 2 emin , ( 13b
)
where • denotes the function that rounds the reals from [0, 2) in the same direction as • but on the regular grid {i • 2 1-p : i = 0, 1, . . . , 2 p }.

In order to handle the cases (13a) and (13b) simultaneously, let us define µ = max(c, e min -2e).

(

) 14
Since • coincides with • on [1, 2), the correctly-rounded value of the exact square is given in both cases by

•(x 2 ) = •(ℓ) • 2 d , (15a)
where

ℓ = m 2 • 2 -µ and d = µ + 2e. ( 15b
)
By construction one has 0 < ℓ < 2 and e min ≤ d ≤ e max . The property below further gives bounds for the range of µ.

Property 5. One has c ≤ µ ≤ p + ǫ with ǫ = p is odd .

Again, these bounds are tight: x = 1 gives µ = 0, while x = ±α ′ gives c = 0, e ′ = e min -p -ǫ, and then µ = p + ǫ.

An explicit formula for •(ℓ). The above analysis has reduced, for precision p, floating-point rounding of the exact result x 2 to fixed-point rounding of the scaled significand ℓ in (15). Using classical rounding formulae (see for example [3, §8.4.3]), we now give an explicit expression for •(ℓ). Writing ℓ = (ℓ 0 .ℓ 1 . . . ℓ p-1 ℓ p . . .) 2 and using ∨ for logical OR, the guard bit and sticky bit of ℓ are, respectively,

g = ℓ p and t = ℓ p+1 ∨ ℓ p+2 ∨ • • • . ( 16
)
Its correctly-rounded value is then given by the formula

•(ℓ) = (ℓ 0 .ℓ 1 . . . ℓ p-1 ) 2 + b • 2 1-p , (17)
where the round bit b satisfies (see [3, pp. 436-437])

b =      g ∧ (ℓ p-1 ∨ t) if • = RN, 0 if • = RD, g ∨ t if • = RU.
(18)

Implementation for the binaryk format

We detail here how to implement, for x generic and the binaryk floating-point format, the computation of r = •(x 2 ) using k-bit integer arithmetic and logic. We assume x is given by its standard encoding into an unsigned k-bit integer X. Since the result r satisfies (15a), it is known (see for example [13, §2.3.1]) that its standard integer encoding R is given by R = D • 2 p-1 + L, where

D = d + e max -1 (19)
and

L = •(ℓ) • 2 p-1 . Using (17) we obtain L = L + b with L = ℓ • 2 p-1 (20)
and b the sticky bit defined in (18), so that eventually

R = D • 2 p-1 + L + b. (21)
Consequently, computing R amounts to deducing D, L, b from X, which we detail now in a parameterized fashion, that is, for the binaryk format. We illustrate our analysis for the binary32 format, and the corresponding C code (for rounding 'to nearest even') appears in Listing 1. In order to give an idea of the ILP exposed by our approach, this C code has been set out in such a way that line i displays only the expressions that can be evaluated in i cycles with the ST231 latencies and assuming unbounded parallelism.

Computing L. From (15b) and (20) it follows that

L = m 2 /2 1-p+µ (22)
and, therefore, we first need to deduce from X an integer encoding of m, say M , as well as the integer µ.

To produce M , recall that m = (1.m 1 . . . m p-1 ) 2 as x is normal. Since p ≤ k a possible choice is to set up m • 2 k-1 , which can be obtained from (3) by shifting and masking:

M = m • 2 k-1 = (X ≪ w) | 2 k-1 .
For the binary32 format, where w = 8 and k = 32, this corresponds to line 2 in Listing 1; on ST231, this takes 2 cycles and, due to the extended immediate value 2 31 = (80000000) 16 , 3 instruction syllables.

To get µ, recall first that x normal implies e = E -e max . Then, recalling that e min = 1 -e max and applying [START_REF] Walters | Combined unsigned and two's complement hybrid squarers[END_REF],

µ = max(c, F ), F = e max + 1 -2E. ( 23
)
To get the boolean value c, it suffices to remark that [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] and

M = m • 2 k-1 imply c = [M > M 0 ], with M 0 = √ 2 • 2 k-1 .
To get the (possibly negative) integer F , first we extract 2E from X in (3) by using the identity

2E = X ≫ (p -2) & (2 w+1 -2), (24)
and then we subtract 2E from the constant e max + 1. For the binary32 format the computation of c and F appears at line 3 of Listing 1, where (b504f 333) 16 is M 0 for k = 32; on ST231 each of c and F takes 3 cycles, so that µ is eventually obtained in 4 cycles.

Let us now see how to deduce L from M and µ. The property below shows that the k most significant bits of the 2k-bit integer M 2 are enough for that purpose.

Property 6. L = ⌊H/2 µ+w-1 ⌋ with H = ⌊M 2 /2 k ⌋.
Given M the mul instruction computes the higher half H of M 2 , which then, by Property 6, simply has to be shifted right by µ + w -1 in order to yield L. For the binary32 format, this appears at lines 5 and 6 of Listing 1. Since here p = 24 is even, we deduce from Property 5 that µ + 7 is at most 31, which thus agrees with the C99 specification of the bitwise shift operator [8, p. 84]. With the ST231 latency constraints, both H and µ + 7 are computed from X in 5 cycles, so that L is obtained in 6 cycles.

Computing b. We focus here on the most interesting case, rounding to nearest even, for which b = g ∧ (ℓ p-1 ∨ t) with g and t as in (16).

Note first that g is the least significant bit of the integer G = ⌊ℓ • 2 p ⌋ = 0≤i≤p ℓ i 2 p-i and, using a proof similar to that of Property 6, we arrive at

g = G mod 2, G = ⌊H/2 µ+w-2 ⌋.
For the binary32 format, the corresponding C code appears at lines 6 and 9 of Listing 1. On ST231, G will have the same latency as L (6 cycles), and we thus get g in 7 cycles. Since ℓ p-1 is the least significant bit of L we have

ℓ p-1 = L mod 2,
so that we are left with computing the sticky bit t. The next result shows how to recover this bit simply by checking that some lower parts of H and X are nonzero, that is, without computing the lower half of the exact square M 2 .

Listing 1. Computing •(x 2 ) for the binary32 format, • = RN, and x generic. For the binary32 format, this subtraction appears at line 10 of Listing 1. Recalling that on ST231 we get F and µ in, respectively, 3 and 4 cycles, we will thus get D in 5 cycles.

Packing the result. From (21) the integer encoding R of the result satisfies R = D • 2 p-1 + L + b and we have just detailed how to get from X the integers D, L, and b. Moreover, assuming the latency model of the ST231, their respective cost has been shown to be of 5, 6, and 9 cycles. This implies a latency of 6 cycles for D •2 p-1 , and using the parenthesization shown at the last line of Listing 1 we eventually get R in 10 cycles. Thus, when • is RN the overall cost is larger than that of the round bit b by only one cycle.

Some simplifications when • is not RN. When the rounding mode • is RD the round bit b in (18) is zero. Consequently, the instructions involving G, T 1 , T 2 , and b can be suppressed and the last line of Listing 1 replaced with: With these new codes the expected latency of R on ST231 drops from 10 to 7 cycles for • = RD, and from 10 to 9 cycles for • = RU.

Detecting and handling special input

We first have to decide whether input x is special or not, that is, to compute from X the value of C spec in [START_REF] Jeannerod | FLIP 1.0: a fast floating-point library for integer processors[END_REF]. By Definition 1 this condition satisfies

C spec = C small ∨ C large ∨ C nan ( 25
)
with 12 } Once special input have been filtered out, it remains to return, for the given rounding mode •, the standard integer encoding R of the associated result r prescribed by ( 8):

C small = |x| < α ′ , C large = |x| ≥ Ω ′ ,
When • is RN. We deduce from ( 6) and ( 8) that for x special, r must be +0 if |x| < α ′ , qNaN if x is NaN, and +∞ otherwise. Implementing this is then straightforward as it suffices to recall from (4) that, on the one hand 0, 2 k-1 -2 p-1 , and 2 k-1 -2 p-2 are standard encodings of +0, +∞, and qNaNs, and that, on the other hand,

C nan = X & (2 k-1 -1) > 2 k-1 -2 p-1 .
For the binary32 format, the computation of C nan is shown at lines 1 and 2 of Listing 2, while lines 6 to 9 display the computation of R. On ST231, lines 6 to 9 will be ifconverted as shown by the following pseudo-code:

Rlarge_or_nan = slct(Cnan,0x7fc00000,0x7f800000) R = slct(Csmall,0,Rlarge_or_nan)
With a latency of 1 cycle for the 'slct' instruction and since C nan costs 2 cycles, we thus get R for x special in 4 cycles.

When • is not RN. For • = RD the only difference with the previous case is when |x| ≥ Ω ′ and, by( 6) and ( 8), we now have r = max(|x|, Ω). The standard integer encoding of Ω is 2 k-1 -2 p-1 -1, which equals (7f 7f f f f f ) 16 for the binary32 format. Consequently, it suffices to replace line 9 of Listing 2 with: else return maxu(absX, 0x7f7fffff); For • = RU, the specification differs from that for • = RN only in the case where |x| < α ′ , for which we have r = min(|x|, α). Since the standard integer encoding of α is 1, an implementation for the binary32 format follows by simply replacing line 6 in Listing 2 by if (Csmall) return minu(absX, 1); On ST231, R still costs 4 cycles as both max(|x|, Ω) and min(|x|, α) have a latency of 2 cycles, like C nan .

Experimental results obtained on ST231

The C codes detailed in Sections 4 and 5 yield a full implementation of squaring, for the binary32 format and each rounding mode. To check correctness we compiled them with gcc (using a C emulation of the mul, max, maxu, and minu operators as given for example in [13, Appendix B]), and compared with the results of multiplication x × x obtained on an Intel R Xeon R workstation. For each rounding mode this exhaustive test took about five minutes.

We also compiled our C codes with the ST200 compiler, in -O3 for the ST231 processor. The remainder of this section details the performances obtained in this context.

Operator performances

Latency and comparison with general multiplication. The latency on ST231 of our binary32 square implementation is shown in the second column of Table 1. Due to if-conversion, it gives for each rounding mode a number of clock cycles independent of the input value x. The values within square brackets indicate the lowest latencies we can theoretically achieve with the ST231 latency constraints and assuming unbounded parallelism; these best latencies follow from our analysis in Sections 4 and 5 and have the form 1 + L with L the best latency for the generic case. This first experiment shows that the latencies achieved in practice are at most 1 cycle form the best possible ones.

For comparison, the third column of Table 1 displays the latencies of the multiply operator x × y available in the FLIP 1.0 library and optimized for the ST231 [START_REF] Jeannerod | FLIP 1.0: a fast floating-point library for integer processors[END_REF]. As shown in the fourth column, our specialization of this multiply operator into a square operator yields a speedup between 1.75 and 2.3, depending on the rounding mode. Instruction-level parallelism. When designing our algorithms in Sections 4 and 5 we strived to expose as much ILP as we could. To evaluate ILP in practice we use instructions-per-cycle (IPC), which is the parallelism really exposed on the target. As shown in Table 2 it is deduced from the assembly code by dividing the number of instructions by the latency. The IPC achieved is close to the highest ILP reachable within the architectural constraints of the machine, demonstrating a very efficient usage of its resources. Comparison with two non-IEEE variants. To study the impact on latency of relaxing the IEEE 754 specification used so far, we have implemented for each rounding mode a finite-math-only variant and a variant without subnormals.

Finite math only. We assume here that input and output are neither infinity nor NaN, and that overflow does not occur. Hence x now satisfies |x| < Ω ′ . On the one hand, this leaves the generic path unchanged, so that the best possible latencies are the same as for our IEEE version. On the other hand, (25) becomes C spec = C small and the C codes of Section 5 can be simplified accordingly. As the third column of Table 3 shows, in practice this simplification has no impact on latency for RD, while it saves 1 cycle for RN and RU.

No subnormals. Here we assume that x is not subnormal, which means, writing λ = 2 emin for the smallest positive normal number, that x is either NaN, zero, or such that λ ≤ |x|. We also assume that if the exact result x 2 lies in the subnormal range (0, λ) then r = +0 for RN and RD, and r = λ for RU. Since x 2 < λ is equivalent to |x| < λ ′ with λ ′ = √ λ, the specification of this non-IEEE variant can thus be deduced from ( 6) and ( 8) simply by replacing α and α ′ with, respectively, λ and λ ′ .

For the special path, this relaxed specification implies C small = [|x| < λ ′ ] and the proof of Property 8 can be adapted to show that we now have

C small = [2E ≤ e max -1].
Thus, it suffices to replace 102 by 126 at line 3 of Listing 2 and, for RU, to replace 1 by 2 23 in the minu operation. These updates clearly have no impact on the latency.

Concerning the generic path, the case (13b) need not be considered anymore, so that µ in ( 14) is now equal to c. Hence the max operation at line 4 of Listing 1 can be removed and we can replace µ by c at line 6. However, as H still has a latency of 5 cycles, this simplification does not shorten the critical path. This means that the best latencies that we can theoretically achieve with this second non-IEEE variant are the same as for our IEEE version. Furthermore, Table 3 shows that this is true in practice as well.

Application examples

We now apply our fast operator to some square-intensive algorithms in order to study its effect on real applications. After giving a theoretical model of the speedup achievable on ST231, we show practical speedups and how they match that model. All experiments have been done on the ST231 cycle-accurate simulator, and while we focus on rounding 'to nearest even' (RN) similar results hold for RU and RD.

Speedup model for loop nests involving squares. While Table 1 gives speedups for a single replacement of multiply by square, we now evaluate the theoretical expectation of speedup for loops. When the square operator ap-pears in a loop of n iterations, we introduce the definition

theoretical speedup = (L + ∆ • σ) • n + C L • n + C
with L, ∆, σ, and C given as follows:

• L is the latency of the loop body where squares are used, which includes the application-related operations inside the loop and the cost to maintain the loop. • ∆ is the latency gap between multiply and square.

• σ is the number of squares in a single iteration.

• C is the cost of the straight-line code outside the loop.

Note that when n tends to infinity, the theoretical speedup tends to 1 + ∆ • σ/L, which does not depend on C.

On ST231 the values ∆, L, and C have the following features. First, Table 1 gives ∆ = 21 -12 = 9 cycles. Second, L and C can be modelled as

L = a • Br + Idx + Ld + C in , C = St + C out .
Here, Br is the cost of a taken branch, which is 3 cycles. The unrolling transformation done by the compiler has the effect of dividing the number of branches taken to jump back to the head of the loop by the unrolling factor; a in (0, 1] is used to denote this effect. The value of a depends on the application and on the compiler optimization level; to estimate the trend of the speedup, we take a = 0.5, meaning that the loop is always unrolled by a factor of 2. Idx is the cost to test loop index, which is 1 cycle; Ld is the cost to load input values, which is 3 cycles. C in (resp. C out ) is the latency of the application-related code within (resp. outside) the loop; for each application, the values of C in and C out can be deduced from the latencies of the 5 basic operators of FLIP 1.0 [13, Table 1] and from the latency of 12 cycles of our square operator. Finally, St is the cost of stack handling of the function call; it ranges from 10 to 12 cycles depending on the achievable ILP of each application.

The above model has been applied to three application examples, which we review now.

Example 1: Euclidean norm. Given a vector v of n floating-point data v i we consider the computation of its Euclidean norm v 2 = ( n i=1 v 2 i ) 1/2 by means of three different algorithms.

Naive algorithm. Here v 2 is produced by a for loop whose ith iteration simply squares v i and adds it to the current partial sum of squares. Figure 1 shows that the theoretical speedup tends to about 1.185 and that in practice we are close to this value as soon as n ≥ 20, assuming n is known at runtime. If n is known at compile time then an even higher speedup is observed, since the compiler achieves more efficient loop unrolling optimization. The nine black bullets in Figure 1 illustrate this fact for n = 2, . . . , 10. The greatest speedups are for n ≤ 4, since in this case all the Two-pass algorithm. This algorithm, which is already mentioned in [START_REF] Blue | A portable Fortran program to find the Euclidean norm of a vector[END_REF], aims to avoid overflow by first computing v ∞ = max i |v i | and then applying the naive algorithm to the scaled vector [v i / v ∞ ] i . The input is scanned twice and n divisions are used. Thus, the speedup we can expect is lower than for the naive algorithm, as shown in Figure 2. However, we see that the practical speedup still matches well the theoretical model as soon as n ≥ 20. One-pass algorithm. This algorithm also intends to avoid overflow but requires only one pass over the data, the scaling factor being now computed on the fly. It is an adaptation of Blue's algorithm [START_REF] Blue | A portable Fortran program to find the Euclidean norm of a vector[END_REF]

. . , v n is 1 n-1 n i=1 (v i -v) 2 , where v = 1 n n i=1 v i .
It is known [6, p. 11] that it can be evaluated accurately by the naive way, which requires two passes over the data: get v first and then iteratively square and add the v i -v's. This method has the same structure as the two-pass algorithm in the previous example (with the maximum replaced by a sum, and the division replaced by a subtraction). Consequently, the theoretical and practical speedups brought by our fast squarer are similar to those in Figure 2.

Example 3: binary powering. For a floating-point datum x and for n an integer power of two such that n ≥ 4, we consider here the evaluation of x n by means of log 2 n successive squares. The results obtained for this application are shown in Figure 3. The qualitative analysis done for the naive algorithm of Example 1 still applies but, since binary powering does not involve any operation other than squaring, higher speedups are observed.

Conclusion and perspectives

This work has presented optimized C codes for IEEE bi-nary32 squaring which are, on ST231, significantly faster than the corresponding implementation of general multiplication, thus demonstrating the practical benefit of operator specialization in this context.

To achieve this, we strived to exploit the specific features of the square operator: it is univariate, always nonnegative, and has predictable overflow ranges; subnormal input can always be handled by the special path; and, in the generic path, the lower half of the exact square of the input significand is not needed, and this even for computing the sticky bit used to provide correct rounding. We also heavily relied on some features of our target processor and its compiler, like VLIW 4-issue parallelism, if-conversion and fast 'slct' instruction, fast min and max instructions, and the availability of extended immediates and of an integer multiplier of type '32 × 32 → 32 high.' However, as our C codes show, the subset of the ST231 instruction set needed by our fast squarer is extremely modest and, unlike for multiplication, the leading-zero count instruction is useless in this context. Furthermore, all our algorithmic descriptions are parameterized by the binaryk format, thus offering increased confidence and flexibility: the correctness analysis is done once for all, and one can deduce portable C codes for various formats in a direct way as soon as the corresponding k-bit integer arithmetic and logic are available. In fact, we already have such a code for the binary64 format and it relies on a 64-bit integer layer highly optimized for ST231. This code is certainly suboptimal, though, as some parts of the algorithm for binary64 squaring could directly use 32-bit integers. Therefore, it would be interesting to implement an optimized 32-bit version of binary64 squaring and to evaluate the latency improvement compared to the 64-bit version produced with our approach.

Three other extensions would also be worth considering. First, although our approach assumes radix 2 we should investigate to which extent it carries over other radices, and especially radix 10. Second, the overhead (in terms of latency and code size) of setting status flags remains to be studied. Third, squaring is not only a specialization of multiplication but also of integer powering, that is, of pown(x,n) = x n with n an integer. Since the 2008 revision of the IEEE 754 standard recommends correct rounding for this operator [7, Table 9.1], its optimized implementation for VLIW integer processors will be a natural extension of the work we have presented here.

A. Proofs

Proof that 2 ≤ p < e max for all standard binary formats. The standard binary formats are defined in [7, Table 3

.5].

There we see that, on one hand, this is true when k ≤ 128. On the other hand, when k > 128 we have p = k -j k + 13 and e max = 2 j k -14 -1, with j k = round(4 log 2 k), the integer closest to 4 log 2 k. (If k is an integer then 4 log 2 k cannot be exactly halfway two consecutive integers, so that no tie can occur.) By definition of round, 4 log 2 k -1/2 < j k < 4 log 2 k + 1/2. This implies first that, for k > 2 7 , 2 ≤ p ≤ k. Second, 2 -14.5 k 4 -1 < e max . Third, k 3 > 2 15.5 and then 2 -14.5 k 4 > 2k > k + 1, so that k < e max .

Proof of Property 1. Since α ′ and Ω ′ in (7) are integer powers of two it suffices to verify that e min ≤ ⌊(e min -p)/2⌋ < (e max + 1)/2 ≤ e max .

The leftmost inequality is equivalent to e min ≤ (e min -p)/2 because e min is an integer; since e min = 1 -e max , the latter inequality is itself equivalent to p < e max , which is true by (2b). From (2b) it also follows in particular that e max ≥ 1, which implies the rightmost inequality. The remaining inequality follows from the fact that (2b) implies that e min is negative while p and e max are positive.

Proof of Property 2. If |x| ≥ Ω ′ then x 2 ≥ 2 emax+1 , so that •(x 2 ) = max • for each •. Conversely, assume that |x| < Ω ′ . Since x is a finite floating-point number in precision p, we deduce that |x| ≤ (2 -2 1-p ) • 2 (emax-1)/2 and then x 2 ≤ C • 2 emax , with C = 2(1 -2 -p ) 2 . Since C < 2 -2 1-p one has further x 2 < Ω. This implies •(x 2 ) < max • for each • and the conclusion follows.

Proof of Property 3. If α ≤ x < α ′ then 0 < x 2 < α/2, so that RN(x 2 ) = RD(x 2 ) = +0 and RU(x 2 ) = α. This shows that x ∈ [α, α ′ ) implies •(x 2 ) = min • for all •. To prove the maximality of α ′ it suffices to check that RN(y 2 ) = min RN for some floating-point number y in [α ′ , 2α ′ ), say y = 3 2 α ′ . We have y 2 = 9 4 2 2⌊(emin-p)/2⌋ and, using the fact that ⌊i/2⌋ ≥ (i -1)/2 when i ∈ Z, we deduce that y 2 ≥ 9 8 α/2 > α/2. It follows that RN(y 2 ) = α, which differs from min RN = +0.

Proof of Property 4. Recalling [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] and taking squares in [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], we obtain 2 2⌊(emin-p)/2⌋ ≤ m ′ • 2 e ′ < 2 emax+1 . On the one hand, m ′ < 2 implies 2⌊(e min -p)/2⌋ < e ′ + 1 and, since both sides are integers, the announced lower bound follows. On the other hand, 1 ≤ m ′ implies e ′ < e max + 1 and, similarly, we deduce the announced upper bound.

Proof of Property 5. The lower bound is an immediate consequence of the definition of µ in [START_REF] Walters | Combined unsigned and two's complement hybrid squarers[END_REF]. To establish the claimed upper bound we consider two cases:

• If µ = c then µ is at most 1 and cannot be larger than p + ǫ, since p ≥ 2 and ǫ ≥ 0.

• If µ = e min -2e then, recalling that e ′ = c + 2e and noticing that the lower bound in Property 4 equals e min -p + ǫ (because e min is even), we obtain µ ≤ p + ǫ + c. Since both µ and p + ǫ are even integers, and since c is either 0 or 1, it follows that µ ≤ p + ǫ. Hence µ ≤ p + ǫ in both cases, and the proof follows.

Proof of Property 6. From (22), M = m • 2 k-1 , and k = w + p it follows that L = ⌊y/n⌋ with y = M 2 /2 k and n = 2 µ+w-1 . Since w ≥ 3 and since, by Property 5, µ ≥ 0, n is a positive integer. To conclude it suffices to apply the fact that ⌊y/n⌋ = ⌊⌊y⌋/n⌋ for n > 0 (see [5, p. 72]).

Proof of Property 7. Let q be the number of trailing zeros of m = (1.m 1 . . . m p-1 ) 2 . Then m 2 can be written m 2 = (s -1 s 0 .s 1 . . . s 2p-2q-2 ) 2 with s 2p-2q-2 = 1.

Thus, H = (s -1 s 0 . . . s k-2 ) 2 and, using (15b) and ( 16), we can also decompose the sticky bit as t = t 1 ∨ t 2 with

t 1 = s p-µ+1 ∨• • •∨s k-2 and t 2 = s k-1 ∨• • •∨s 2p-2q-2 .
By Property 5 and since w ≥ 3, we have k -2 -(p -µ + 1) + 1 = µ + w -2 ∈ {1, . . . , k -1}. Hence t 1 = 1 if and only if the last µ + w -2 bits of H are not all zero, that is, if and only if the integer T 1 obtained by shifting H left by k -(µ + w -2) = p + 2 -µ is nonzero. Since s 2p-2q-2 = 1 we have t 2 = 0 if and only if k -1 > 2p -2q -2, that is, if and only if the number q of trailing zeros of m is at least p -⌊k/2⌋. The latter condition is equivalent to X mod 2 p-⌊k/2⌋ = 0, which concludes the proof.

Proof of Property 8. Let |X| and A ′ denote the standard integer encodings of |x| and α ′ , respectively. It is known [10, p. 58] that |x| < α ′ if and only if |X| < A ′ , that is, |X| ≤ A ′ -1.

(26) By ( 3), |X| = (E + ǫ) • 2 p-1 with ǫ ∈ [0, 1 -2 1-p ] and, by Property 1, A ′ = E ′ • 2 p-1 with E ′ = ⌊(e min -p)/2⌋ + e max . Thus, (26) is equivalent to E ≤ E ′ -(ǫ + 2 1-p ). Since E, E ′ are integers and ǫ + 2 1-p ∈ (0, 1], the latter inequality is equivalent to E ≤ E ′ -1. Now, e min = 1 -e max gives E ′ -1 = ⌊(e max -p -1)/2⌋ and we conclude using the fact that i ≤ ⌊j/2⌋ is equivalent, for integers i, j, to 2i ≤ j.

Proof of Property 9. We use the same notation as in the proof of Property

Property 4 .

 4 The normalized exponent e ′ of x 2 satisfies 2⌊(e min -p)/2⌋ ≤ e ′ ≤ e max .

1Property 7 .

 7 T2 = X & 0xff; 2 M = (X << 8) | 0x80000000; E2 = (X >> 22) & 0x1fe; 3 F = 128 -E2; c = M > 0xb504f333; 4 mu = max(c, F); 5 H = mul(M, M); 6 L = H >> (mu + 7); G = H >> (mu + 6); T1 = H << (26 -mu); 7 8 9 b = (G & 1) && ((L & 1) | (T1 | T2)); 10 return (((mu -F) << 23) + L) + b; One has t = [T 1 = 0] ∨ [T 2 = 0] with T 1 and T 2 the two k-bit integers given byT 1 = H ≪ (p + 2 -µ) and T 2 = X mod 2 p-⌊k/2⌋ .For the binary32 format, Property 7 gives T 2 = X mod 2 8 , which can be implemented by masking X as shown at line 1 of Listing 1. The computation of T 1 is a mere left shift of H by 26 -µ, the latter value ranging in [0, 31] thanks to Property 5. Then notice that the bit ℓ p-1 ∨ t is zero if and only if the integer U obtained by bitwise-ORing the integers L & 1 and T 1 andT 2 is zero. The logical AND of g = G & 1 ∈ {0, 1} and U is thus enough to yield b, which allows us to avoid testing explicitly if T 1 or T 2 is nonzero. This is shown at line 9 of Listing 1. The parenthesization chosen there aims to reduce the overall latency for b on ST231: both L and T 1 can be obtained in 6 cycles, while T 2 costs 1 cycle; therefore, both L & 1 and T 1 | T 2 follow in 7 cycles, which yields b in 9 cycles. Computing D. From the definitions of d and D in (15b) and (19) we deduce that D = µ + 2e + e max -1. Hence, recalling that e = E -e max and the definition of F in (23), D = µ -F.

  return ((mu -F) << 23) + L; When • is RU the bit ℓ p-1 is not needed and b is the logical OR of g = G & 1 and T 1 | T 2. In this case, we thus replace line 9 of Listing 1 by: b = (G & 1) | (T1 | T2);

Figure 1 .

 1 Figure 1. Impact of square on naive Euclidean norm.

Figure 2 .

 2 Figure 2. Impact of square on two-pass Euclidean norm.

  attributed to Hammarling and implemented in LAPACK [6, p. 507]. Each of the n iterations performs exactly 1 square as well as 2 or 4 additional operations, depending on the data v i . The total number of operations thus varies dynamically and turns out to be maximum when |v 1 | < |v 2 | < • • • < |v n |, and minimum when |v 1 | ≥ |v 2 | ≥ • • • ≥ |v n |. Each of these two extreme cases yields a behavior similar to the one displayed in Figure 2, the limiting value when n → ∞ being about 1.077 in the first case, and about 1.096 in the second case. Example 2: sample variance. The sample variance of v 1 , .

Figure 3 .

 3 Figure 3. Impact of square on binary powering.

  For the binary32 format, a C fragment implementing C spec by means of 2E and the two previous properties is shown at lines 2 to 4 of Listing 2. On ST231 the cost will be of 4 cycles and, as C spec is independent of •, this fragment holds not only for • = RN but also for • ∈ {RD, RU}.

	Listing 2. Detecting and handling special in-
	put for the binary32 format and • = RN.
		absX = X & 0x7fffffff;	1
	E2 = (X >> 22) & 0x1fe; Cnan = absX > 0x7f800000;	2
	Csmall = E2 <= 102;	Clarge_or_nan = E2 >= 382;	3
	Cspec = Csmall || Clarge_or_nan;		4
	if (Cspec) {			5
	if (Csmall) return 0;	// r = +0	6
	else {			7
	if (Cnan) return 0x7fc00000;	// r = qNaN	8
	else return 0x7f800000; }	// r = +oo	9
	} else {			10
	// generic case (Listing 1).		11

and C nan = x is NaN . The next two properties show how to obtain C small and C large ∨ C nan by reusing the value 2E computed for the generic case (see (24) and line 2 of Listing 1). Property 8. C small = [2E ≤ e max -p -1]. Property 9. C large ∨ C nan = [2E ≥ 3e max + 1].

Table 1 .

 1 Latency comparison for square and multiply.

	•	square	FLIP 1.0 multiply	speedup
	RN	12 [11]	21	1.75
	RD	9 [8]	21	2.3
	RU	11 [10]	21	1.9
	RZ	9 [8]	18	2

Table 2 .

 2 Latency, code size, and IPC for square.

	•	Latency L	Number N of instructions	IPC = N/L
	RN	12	42	3.5
	RD	9	31	3.4
	RU	11	37	3.4

Table 3 .

 3 Latency comparison with two non-IEEE variants.

  8 and write O ′ for the standard integer encoding of Ω ′ . The special integer encoding used for NaNs gives C large ∨ C nan = [|X| ≥ O ′ ]. By Property 1 we haveO ′ = (3e max + 1)/2 • 2 p-1 , so that |X| ≥ O ′ is equivalent to E + ǫ ≥ (3e max + 1)/2. Since ǫ ∈ [0,[START_REF] Bertin | Techniques and tools for implementing IEEE 754 floating-point arithmetic on VLIW integer processors[END_REF], this is equivalent to E ≥ (3e max + 1)/2 and the conclusion follows.