Aymeric Hervieu
email: aymeric.hervieu@inria.fr

Benoit Baudry
email: benoit.baudry@inria.fr

Arnaud Gotlieb

PACOGEN : Automatic Generation of Pairwise Test Configurations from Feature Models

Feature models are commonly used to specify variability in software product lines. Several tools support feature models for variability management at different steps in the development process. However, tool support for test configuration generation is currently limited. This test generation task consists in systematically selecting a set of configurations that represent a relevant sample of the variability space and that can be used to test the product line. In this paper we propose PACOGEN to analyze feature models and automatically generate a set of configurations that cover all pairwise interactions between features. PACOGEN relies on constraint programming to generate configurations that satisfy all constraints imposed by the feature model and to minimize the set of the tests configurations. This work also proposes an extensive experiment, based on the state-of-the art SPLOT feature models repository, showing that PACOGEN scales over variability spaces with millions of configurations and covers pairwise with less configurations than other available tools.

I. INTRODUCTION

Feature models (FMs) allow companies to reason over a large number of variants for their software systems [START_REF] Benavides | Automated analysis of feature models: A detailed literature review[END_REF]. Several tools exist to manage variability with feature models [START_REF] Beuche | Variability management with feature models[END_REF], debug configurations [START_REF] White | Automated diagnosis of feature model configurations[END_REF] and to derive specific configurations [START_REF] Czarnecki | Staged configuration using feature models[END_REF]. However, very few techniques or tool support the testing activity in software product line engineering. The size of the variability space is a major challenge for testing. Realistic feature models can represent millions of variants, which means that exhaustive testing is impossible in most cases. For example, the feature model for Arcade video games in the SPLOT repository [START_REF] Mendonc ¸a | [END_REF] is composed of 61 features, but more than 1 million valid products can be derived from this model.

One challenge consists in selecting a small subset of all possible configurations for testing. This sample should cover relevant characteristics of the feature model, while staying as small as possible. Keeping a small sample is crucial to limit the effort necessary for testing each selected configuration. Recent work suggest exploring pairwise coverage [START_REF] Cohen | The aetg system: An approach to testing based on combinatorial design[END_REF] to sample the variability space. The automatic generation of a minimal set of configurations that cover pairwise interactions faces two challenges: (i) feature models specify dependencies between features which forbid some pairwise feature interactions; (ii) the generation of a minimal pairwise test set is a complex optimization problem. Perrouin et al. [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] and Oster et al. [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] have recently investigated the first issue. In this work, we explore both issues at the same time (minimization and dependencies between features).

We propose PACOGEN 1 for the automatic generation of test configurations that cover all valid pairwise interactions in a feature model. PACOGEN has two major characteristics: it processes feature models directly from Eclipse, the most common IDE for model-driven development; test generation is based on constraint programming. We choose constraint programming first because of its flexibility to deal with dependencies between features for test generation: this programming paradigm allows us to design a tailor-made pairwise constraint. Second, unlike SAT-solving, constraint programming is well suited for optimization problems such as those related to the minimization of the size of test sets. PACOGEN users can decide to ask for the smallest test set that covers pairwise interactions, or they can ask for the smallest solution that can be found in a given amount of time (anytime minimization).

The paper presents the following contributions: 1) PACOGEN , a constraint-based testing tool for automatic generation of test configurations that cover all pairwise interactions in a feature model ; 2) A series of experiment with 69 feature models from SPLOT 2 , one of the largest and up-to-date feature model repository. The main results show that our strategy for pairwise generation scales over variability spaces that specify millions of configurations and that we can generate less configurations than state of the art techniques [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF], [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF]. The rest of the paper is organized as follows. Section 2 presents some background on combinatorial interaction testing. Section 3 details our constraint-based model of FMs and how pairwise coverage can be enforced in a set of test configurations. Section 4 explains how the number of test configurations can be minimized through the usage of well-known Constraint Programming techniques. Section 5 presents experiments that compare our approach to other techniques and that run the generation 69 models from SPLOT. Finally, Sec. 6 draws some conclusions and perspectives of this work.

II. BACKGROUND

This section briefly introduces the metamodel we use to build feature models. Then, we define pairwise interaction coverage over a feature model. We also emphasize some issues that must for the automatic generation of a minimal test configurations set.

A. Feature models

Perrouin et al. [START_REF] Perrouin | Reconciling automation and flexibility in product derivation[END_REF] have built a metamodel that formally captures the definition of a feature model, on the basis of the work by Schobbens et al. [START_REF] Schobbens | Generic semantics of feature diagrams[END_REF]. The metamodel is displayed in figure 1 and defines the structure of a feature model as follows:

• A feature model (FEATUREDIAGRAM class) is composed of a set of FEATUREs. • We distinguish between FEATUREs and PRIMITIVE-FEATUREs. Our test generation process considers only interactions between primitive features. • One parent FEATURE is related to a set of children features through an OPERATOR or through a binary constraint. A FEATURE can also have a list of AT-TRIBUTEs.

• Five different OPERATORs can relate a parent feature to its children: AND, OR, XOR, OPT, CARD. • Two CONSTRAINTs can relate features that are not parent / children: REQUIRE and MUTEX. Figure 2 displays a small feature model for a car break system (White et al. [START_REF] White | Automated reasoning for multi-step feature model configuration problems[END_REF]). This model conforms to the metamodel definition. It specifies that a car has a backward sensor that can associated with a Lateral Range Finder (LRF) or Forward Range Finder (FRF). A car also has an optional Automated Driving Controller (ADC). An ADC must have a Collision Avoidance Breaking (CAB): either Standard Avoidance (SA) or Enhanced Avoidance (EA). An ADC also has an option for Parallel Parking (PP). The feature model specifies two cross-tree constraints: if the PP option is chosen, then the car must have a LRF; if the EA is chosen, then the car must have FRF.

B. Pairwise testing for selecting test configurations

Pairwise testing is a particular case of combinatorial interaction testing (CIT) introduced by Cohen et al. [START_REF] Cohen | The aetg system: An approach to testing based on combinatorial design[END_REF] to sample large test input domains. This test selection technique focuses on the subset of the input domain that covers all value combinations for each pair of variables. For example, the car configuration [ADC, EA, PP, LRF, FRF] covers the interactions between the following pairs of values: For pairwise testing over feature models, the strength t of the array is equal to 2, v is the number of primitive features and all these variables are defined on a domain v i of size 2 (all boolean variables).

A general issue for pairwise test generation is to determine the number N of necessary lines in the array in order to cover all interactions between variable pairs. This is a complex optimization problem that has been tackled mainly through three types of solutions [START_REF] Cohen | Interaction testing of highly-configurable systems in the presence of constraints[END_REF]: algebraic constructions, greedy algorithms and meta-heuristics. However, most previous solutions consider that all values for all variables are independent, i.e., all value interactions are possible.

Building a set of pairwise configurations from feature models requires to think about new generation strategies because a number of value interactions are forbidden. For example, the pair (EA = selected, FRF = unselected) is not a valid pair according to the feature model in figure 2. Thus the optimization problem we tackle here can be formulated as follows

Given a feature model we look for the minimum number of valid configurations that cover all the authorized interactions between pairs of features. In the following section we introduce a new constraint model for feature model that targets both the generation of valid configurations with respect to the feature model and the minimization of test configurations to cover pairwise interactions.

III. FEATURE MODELS AS CONSTRAINTS

Our approach is based on constraint programming through the mapping of a feature model into a finite domain constraint model. This section details our model and the structure used to generate test configurations for pairwise testing. We start by recalling what is a global constraint (Sec.III-A). Then, we detail the constraints generated for modeling parent/children operators from the feature model and those generated for handling cross-tree links (Sec.III-B). Finally, we explain how to enforce pairwise coverage on a set of test configurations (Sec.III-C). An important contribution of our Feature Model of a Car Break System [START_REF] White | Automated reasoning for multi-step feature model configuration problems[END_REF] work is the definition of a new global constraint that enforces a set of values to be included within a pair of variable vectors. This constraint was not reported elsewhere, although the usage of global constraints for handling feature models was already skeched [START_REF] Karatas | Global constraints on feature models[END_REF]. The filtering algorithm developed for this global constraint is a key point of our approach as it allows pairwise coverage to be enforced within a set of test configurations.

A. Global constraints

An important feature of constraint programming techniques is their ability to allow users to define new specialpurpose constraints under the form of global constraints. A global constraint is a relation defined by an interface (operator's name and constrained variables), a filtering algorithm and awakening conditions. For example, the global constraint all different([X 1 , .., X n ]) constrains the set of finite domain variables X 1 , .., X n to take distinct values whatever those values. Its filtering algorithm is based on the matching theory [START_REF] Régin | A filtering algorithm for constraints of difference in csps[END_REF] and its awakening conditions may vary from one implementation to another but often correspond to the discovery that one of the variables has become instantiated. In fact, a global constraint is handled as any primitive constraint by the solver and its filtering algorithm is launched anytime one of its awakening conditions is satisfied during constraint propagation.

B. Constraints for feature model links

From a feature model, the problem of generating a sequence of test configurations enforcing pairwise coverage is mainly a constraint modelling problem. The model includes two sets of constraints: the constraints handling inheritance links and the constraints handling cross-tree links. As described above, inheritance links are hierarchical relations between features (namely, OPT, AND, OR, XOR, CARD) while cross-tree links are alternative relations among unrelated features (namely, REQUIRE, MUTEX). Based on their definition, we have defined dedicated global constraints to capture these relations. For inheritance links, global constraints

opt(A, [B, C, ...]), and(A, [B, C, ...]), or(A, [B, C, ...]) xor(A, [B, C, ...], card(A, [B, C, ...], N, M)
holding over two-valued (i.e., 0 and 1) domain variables are used. These relations express the hierarchical relation among the feature parent A and his children B, C, .... The cardinality global constraint existed already in the CP community [START_REF] Hentenryck | The cardinality operator: A new logical connective for constraint logic programming[END_REF]. In addition to the A, B, C, ... variables, it takes two additional variables N and M as inputs for the minimum and maximum number of children, but these variables are always instantiated in feature models. For REQUIRE and MUTEX, we encoded the corresponding logical relations (implication and mutual exclusion respectively) onto the two-valued domain variables associated to feature children. No further details on these relations are necessary as they are trivial to implement.

C. Enforcing pairwise coverage through constraints

In order to build a sequence of test configurations enforcing pairwise coverage of a given feature model, we defined a specific data structure that is incrementally filled in. The data structure is basically a matrix with columns containing values of all the features and rows representing test configurations. Our framework fills in this matrix in order to find a sequence of test configurations (the rows of the matrix) that covers the pairwise criterion. At the begining of the process, the matrix is filled in with two-valued unknown variables and its size is positionned to a sufficently large constant 3 as the number of test configurations is unknown (and cannot be computed from analytical results [START_REF] Hnich | Constraint models for the covering test problem[END_REF]). At the end of the process, this matrix encodes a solution of the problem, that looks like:

Conf. F eat. A B C ... 1 1 0 0 ... 2 0
I becomes instantiated to 2 then (X 2 , Y 2 ) = (1, 1)
whereas if X 3 is instantiated to 0 then 3 will be removed from the domain of I. In this latter case, the pair (X 3 , Y 3 ) cannot be equal to (1, 1) but there is not enough information to instantiate the variables. The constraint pairwise will be suspended until more information becomes available. The filtering algorithm shown below is used for constraint pairwise(I, (L 1 , L 2 ), (v 1 , v 2 )), where I is finite domain variable, L 1 and L 2 are two lists of finite domain variables, and (v 1 , v 2 ) is a pair of values. This algorithm is launched 3 In practice, we used the value 50 Input: I a finite domain, L 1 , L 2 two lists of finite domains of the same size and (v 1 , v 2 ) a pair of integer values Output: F ail or pruned domains for (I, L 1 , L 2 )

function pairwise(I, (L 1 , L 2 ), (v 1 , v 2 )) I ← I, T 1 ← ∅, T 2 ← ∅; foreach i ∈ I do if (v 1 ∈ L 1 [i]) or (v 2 ∈ L 2 [i]) then I = I \{i} else T 1 ← T 1 ∪ L 1 [i], L 1 [i] ← L 1 [i]; T 2 ← T 2 ∪ L 2 [i], L 2 [i] ← L 2 [i]; if I = {a} then L 1 [a] = v 1 ; L 2 [a] = v 2 ; return ({a}, L 1 , L 2 ); else if (I = ∅ or v 1 ∈ T 1 or v 2 ∈ T 2 ) then return Fail else return (I , L 1 , L 2 )
each time at least one of the domains of

I,L 1 [i] or L 2 [j]
is pruned (awakening conditions). The underlying idea is to explore each of the possible values of I and to determine whether this value is still consistent with the domain of other variables. The complexity of this algorithm is linear w.r.t. the domain size of I, as it iterates only on the possible values of I.

IV. PAIRWISE COVERAGE AS A TIME-CONSTRAINED MINIMIZATION PROBLEM

As said previously, finding the minimum number of test configurations covering the pairwise criterion is a challenging problem. This section explains how this problem can be addressed using the constraint model described above. We explain first how to formulate it as an optimization problem (Sec.IV-A), and second, we solve it using an anytime labelling search procedure (Sec.IV-B).

A. An optimization problem

Our goal is to find the minimum number of test configurations, i.e., the minimum number of lines to instantiate in the matrix for covering the pairwise criterion. This can be achieved by searching the minimum of a cost function f , as follows: Find I 1 , ..., I 4n 2 such that Min(f ) And ∀i, j in 1. Several functions can be considered for minimization. In our framework, we explored two semantically equivalent formulations, namely f 1 = k∈1..4n 2 I k and f 2 = Max k∈1..4n 2 I k Both functions can be used in order to find the minimum number of values for the I k , such that pairwise is satisfied in the matrix. For solving the optimization problem, we used the well-known branch-and-bound method that explores feasible solutions while maintaining the cost function as low as possible. Roughly speaking, at each node of the search tree, the branch-and-bound method evaluates the cost function, prunes subtrees for which the cost will be clearly higher than a current value and selects the subtree that has the least cost. Several parameters impact the search, including the way variables and values are selected for labelling. Another characteristic of the constraint solving techniques we used is their versatility for addressing feature models. This has already been reported in the literature [START_REF] Salinesi | Exploiting the versatility of constraint programming over finite domains to integrate product line models[END_REF]. In particular, additional constraints can easily be defined to take into account extended feature models and several search heuristics can be exploited to generate test configurations. In our framework, we selected the variable with the smallest domain to be enumerated first and from the domain of this variable, we selected the smallest value first. Other heuristics may be chosen but our experimental results showed that these ones are sufficient to handle the largest feature models.

B. Anytime minimization

The constraint solving techniques used in our framework share an interesting property with anytime algorithms [START_REF] Zilberstein | Using anytime algorithms in intelligent systems[END_REF]: they can be stopped at any time or can be given a time or resource contract. In our framework, we gave a time-contract to the branch-and-bound method. As a result, as soon as a first feasible solution has been found, it returns the optimal number of test configurations enforcing pairwise coverage, found only in a given amount of time (ranging from a few seconds to more than three hours). Of course, better feasible solutions might be found if more time is allocated to the search. However, allocating a time-contract to the search permits one to balance advantageously between quality of the solution w.r.t. time needed to find it.

In our experiments, we observed that the branch-andbound method computes good-quality solutions (i.e., nearoptimal solution) in little amount of time and most of the remaining time is used to prove that no better quality solutions really exist. Hence, by relaxing the problem to the finding of near-optimal solutions only, we got a very efficient way to solve our challenging minimization problem.

V. PACOGEN

In this section, we present our PACOGEN implementation and we show how it processes a feature model for generating a set of valid pairwise-covering test configurations (i.e., a set of configurations that covers all pairwise interactions between features). PACOGEN has characteristics that make it suitable for software product line engineering (SPLE). First, the input Feature Model is interpreted as an instance of the metamodel of Fig. 1. This design decision is meant to include PACOGEN directly in a modelling environment dedicated to SPLE. In addition, the processing of Feature Models is independant from the test configurations generation process Second, we used constraint programming to develop PACO-GEN , which is flexible enough to easily customize the test configurations generation process.

Implementation. The test configurations generation process shown in Fig. 3 can be explained as follows: first, PACOGEN transforms the feature model into a constraint model ; second, it adds all the valid pairs under the form of pairwise constraints to the constraint model and then it fills in a special data structure, called constrained matrix, with constrained variables ; third, all the variables in the matrix are labelled in order to satisfy pairwise coverage. We implemented this process within the PACOGEN tool that is freely available 4 . PACOGEN contains four main components:

1) FM Analyzer: this component transforms a feature model into a constraint model, under the form of an abstract syntax tree. As an example of the concrete structure generated within PACOGEN , consider the constraint model generated for the Car Break System FM example of Fig. 2 (called CarFM in the following) that is shown below. FList corresponds to the features, while CList represents the relations, and solver corresponds to the call to the constraint solver with a list of parameters. Size corresponds to an over-estimation of the size of the matrix. This parameter can be used both for proving that there is no solution under a given threeshold or to refine an existing bound. A default value of 50 is usually a good threeshold to start with, but note that the minimization process will always try to find the smallest size within where

CCAR = ⎛ ⎜ ⎜ ⎜ ⎝ CAR 1 CAR 2 . . . CARn ⎞ ⎟ ⎟ ⎟ ⎠ CADC = ⎛ ⎜ ⎜ ⎜ ⎝ ADC 1 ADC 2 . . . ADCn ⎞ ⎟ ⎟ ⎟ ⎠ CCAB = ...

4)

Pairwise test configurations generation: this component calls the constraint solver and the anytime minimization process. It generates the first-found solution under the form of an instantiated constrained matrix, such as the one shown in Fig. 4. Tool validation. In order to validate our tool PACOGEN , we implemented some automated analysis and checked our tool results with published results:

• we computed the total number of valid test configurations with PACOGEN and checked the results with those provided by the SPLOT repository [START_REF] Mendonc ¸a | [END_REF]. On all the cases but one, PACOGEN gave us the published results. After investigation, it turned out that the FM for which we found a difference was due to an interpretation difference within the cross-tree constraints of SPLOT. We corrected our constraint model to mimic the semantics of SPLOT cross-tree constraints ; • on every SPLOT FM, we checked that all pairwise interactions were actually covered by 1) generating all possible pairs of features without cross-tree constraints, and 2) checking that all pairs uncovered by PACOGEN were indeed invalid pairs (i.e., pairs that do not satisfy the cross-tree constraints specified by the FM). The results showed that PACOGEN behave as expected.

VI. EXPERIMENTAL EVALUATION

This section introduces a series of experiments to evaluate Pacogen w.r.t. available techniques and tools. We performed these experiments on the SPLOT repository [START_REF] Mendonc ¸a | [END_REF], which contains more than a hundred FMs. SPLOT has been developed to support empirical studies on the performance and scalability of automated techniques for reasoning on FMs. All the experiments were performed on a standard Intel Core i7 CPU 2.67GHz with 4GB memory.

A. Experimental results

The goal of our first experiment was to evaluate the capability of PACOGEN to find the minimum set of test configurations that covers all the pairwise interactions among features. Although the approach of Oster et al. [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] work section), its greedy algorithm also aims at minimizing the number of configurations and it is currently considered as the most efficient approach for this task [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF]. Tab.I shows the experimental results we got with PACOGEN and how it compares with the results of [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] on ten FMs from SPLOT. Tab.I contains the number of features of each FM (#F), the number of valid configurations (#N) to characterize the size of the search space, the results of both [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] and PACOGEN in terms of number of configurations to cover all the valid pairwise interactions between features, and finally the gain obtained with PACOGEN . Negative gains indicate that PACOGEN proposes less configurations than other approaches, which was the main objective of our work. The goal of this experiment was to evaluate the minimization process and thus we did not try to optimize on the test configuration generation time. In all the cases but one, the generation time required by PACOGEN was less than a few minutes. On Electronic shopping, PACOGEN took 12 hours to get an optimal result which remains an acceptable amount of time for a one-shoot generation. The results indicate that, in average, our approach requires 28.2% lesser configurations than the [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] results. In real-world applications, testing a configuration may be costly because it may require to build the configuration with physical components, to set up a dedicated testing environment and to execute a timeconsuming result evaluation process. Depending on the cost and time required to test a configuration, the improvement we got with PACOGEN may be crucial. Still, there is an anormal result for Inventory where PACOGEN selected more configurations of the results given in [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF]. We manually double-checked the results of PACOGEN without finding how to get a smaller value than 15 for covering all the valid pairwise interactions on this example. Apart from this example, our results show that the PACOGEN approach is well-suited to find the minimum number of configurations that cover all the pairwise interactions of features in a FM.

The goal of our second experiment was to evaluate the anytime minimization process of PACOGEN , which permits to find the best compromise between generation time and result quality. We conducted a large-scale experiment over 67 FMs extracted from SPLOT that is reported in Fig. 5. We launched PACOGEN to compute the minimum number of test configurations by allocating 2 seconds, 10 seconds, 3 min and 15 min to the anytime minimization process. First, the curve shows that PACOGEN can provide an optimal solution (a minimum set of test configurations covering all pairwise feature interactions) in less than 15min in all the cases. Second, the anytime minimization process is interesting to refine the set of configurations in a number of cases. For almost half of the FMs, the anytime minisation process has been useful. Third, the number of test configurations that is refined during this process is usally very small ; there are ony two cases where the number of configurations is refined of more than 5 configurations. These results show that PACOGEN implements a usefull anytime minimization process that has potential to help finding good compromises between number of configurations and generation time.

B. Threats to validity

External threats to validity lie on the source of the empirical data. We have selected our subjects in the SPLOT repository, which allowed us to experiment over a large set of models and to compare our results to Oster's. However, SPLOT is fed with models that mostly come from academia and might not perfectly reflect industrial usage of feature models in terms of number of features and size of the variability space. The threats to internal validity come from biases that can be introduced by errors in the prototypes we used to perform our experiments. In order to limit these threats, we performed exhaustive unit testing on the transformations used PACOGEN in addition to the tests mentionned in Sec. V.

VII. RELATED WORK

Benavides et al. [START_REF] Benavides | Automated analysis of feature models: A detailed literature review[END_REF] surveyed the various constraint models that have been proposed in the literature. In particular, Constraint Satisfaction Problems (CSPs) have been used to encode automated analyses such as filtering the configurations or checking the consistency of the feature model (e.g., dead feature, valid product). Recently, Karatas et al. proposed the usage of global constraints for handling with extended feature models [START_REF] Karatas | Mapping extended feature models to constraint logic programming over finite domains[END_REF], [START_REF] Karatas | Global constraints on feature models[END_REF]. Our approach shares similarities with these works, as our model is also based on finite domain constraints filtering. However, it differs on its goals, our approach being targeted to generate pairwise covering set of configurations, and on its underlying technology. A key-point of our approach is the development of a new global constraint (i.e., pairwise) and the usage of a time-contract optimization labelling procedure, that both are unavailable elsewhere. Another trend in Software Product Line Testing is related to the usage of SAT-solving to generate test configurations [START_REF] Benavides | Automated analysis of feature models: A detailed literature review[END_REF]. The recent work of Uzuncaova et al. [START_REF] Uzuncaova | Incremental test generation for software product lines[END_REF] showed that an extended SAT solver could be used to incrementally generate test inputs for each product in a product line. However, one advantage of Constraint Programming over SAT-solving is its capability to solve optimization problems, instead of pure satisfiability problems [START_REF] Nieuwenhuis | On sat modulo theories and optimization problems[END_REF]. The problem we attacked in this paper is clearly an optimization problem (i.e., to generate the smallest set of configurations ensuring the coverage of the pairwise criterion over the possible feature values).

The main related work that share our testing objectives are the proposals by Perrouin [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] and by Oster [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF]. Perrouin et al. transform feature models in Alloy to select configurations that are valid with respect to the initial model. This approach was the first work to propose the adaptation of pairwise testing in the context of feature modelling, but it faces scalability issues. First because it uses a generate-and-test approach to select the valid pairs that must be covered (generate all possible pairs, then filter all invalid pairs). Second, the Alloy model needs to be transformed in CNF formulae before it can be solved by SAT solvers. This additional step increases the scalability problems. Our approach limits these problems by adopting a test-and-generate approach for identifying valid pairs and by transforming the feature model directly into constraint problem over finite domains. Oster et al. deal with dependencies between features by a flattening transformation over the feature model and the adaptation of the AETG [START_REF] Cohen | The aetg system: An approach to testing based on combinatorial design[END_REF] algorithm in order to consider dependencies between features. The authors do not provide a detailed description of their transformation process, but their experimental results are good. Their approach is efficient and effective: it can select a small number of configurations over large feature models in very reasonable times. These work can be seen as an extension of those of Cohen et al. [START_REF] Cohen | Constructing interaction test suites for highly-configurable systems in the presence of constraints: A greedy approach[END_REF] where a SAT-solver is used to build a greedy approach tackling the problem of dependencies within highlyconfigurable systems. However, nor this approach or Olster's or Perrouin's work tackle the selection of the minimum number of configurations to cover pairwise interactions.

In our work, we propose a deterministic (as opposed to "greedy") approach able to compute the optimal number of configurations. If necessary, our approach can be relaxed by using a time-contract process to provide a near-optimal value. In addition, our work is the first to perform a large experiment on a significant set of feature models .

VIII. CONCLUSIONS AND PERSPECTIVES

In this work we proposed PACOGEN a tool for generating a set of test configurations from FM that covers all the valid pairwise interaction among features. The tool offers an anytime minimization process that allows the user to define an objective stipulating the amount of time he allocates to the generation. Based on Constraint Programming techniques including a dedicated branch-and-bound algorithm, the tool returns the minimum number of configurationd found in the allocated time. We conducted a large-scale experiment over 67 FMs extracted from the SPLOT repository [START_REF] Mendonc ¸a | [END_REF] showing that 1) PACOGEN overcomes the State-Of-the-Art technique of [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] in terms of number of configurations and 2) the anytime minimization process is useful in practice. More precisely, PACOGEN generated 29.2% less configurations than [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] in average and generated the optimal number of configurations for the 67 FMs in less than 15 minutes.

On the basis of these promising results we plan first to continue investigating constraint-based exploration of variability spaces in future work. In particular, we plan to include the values of attributes associated to features in the selection process. For example, if features are associated to execution times, or energy consumption, it is possible to reuse our constraint model to optimize quality of service requirements. Second, we will investigate extensions of our work to consider the relationships between several variability spaces. For example, in component-based systems, the configuration of one component can have an impact on the way another component can be configured. Thus, when sampling the configuration space for the first component it is necessary to consider the dependencies with the configuration space of the second component. This should prevent selecting only configurations that force the second component to be in a narrow area of its configuration space.

THANKS

We are grateful to Pierre Flener and Joseph Scott for some valuable comments and suggestions on an earlier version of this article.

Figure 1 Figure 2 .

 12 Figure 1. A metamodel for Feature Models

  .n, pairwise(I k , (C A , C B ), (1, 1)), pairwise(I k+1 , (C A , C B ), (1, 0)), pairwise(I k+2 , (C A , C B ), (0, 1)), pairwise(I k+3 , (C A , C B ), (0, 0)) where n denotes the number of features while the C A , C B denote the columns of the matrix representing features A and B. Each I k denotes the line of the matrix. Note that additional constraints implicitly enforce all the C A [i], C B [j] to be part of the feature model.

Figure 3 .

 3 Figure 3. The test configurations generation process of PACOGEN

Figure 5 .

 5 Figure 5. Number of test configurations using various time-allocations

  1 and the parameter value. TimeOut and Minimization allow the user to customize the anytime minimization step with its own values. Default values can be used for those who have not any requirement on the test configurations generation time or the size of the test set. 2) Consistency checker: this component evaluates the FM constraints within the constrained matrix data

	Meta
	Model	Pacogen
			FM
			Analyzer
	Conforms to
	Feature
	Model .xml	Consistency checker	Constrained Matrix	Pairwise Constraint generator
		User parameters:
		-Matrix size
		-time_out
		-minimize(rank)	Test	Pairwise(I1,(CCAR,CADC),(0,1))
		minimize(size) -….	configurations generator	Pairwise(I2,(CCAR,CADC),(0,0)) …
		TC1 = (ADC,Car, PP, …)
		TC2 = (Car, PP, …)
		….

  do not explicitely target this objective (see below in the relatedCAR ADC CAB SA EA PP BACKSENSOR LRF FRF

	Figure 4. A constrained matrix solution covering all pairwise interactions of the CarFM
	Feature models	#F	#N	Oster et al. 2010	PACOGEN	Gain
	crisis management	17	-	15	5	-66.7%
	Smart Home	35	1 048 576	11	8	-27.3%
	Inventory	37	2 028 096	12	15	+20%
	Sienna	38	2 520	24	20	-16.7%
	Web portal	43	2 120 800	26	16	-38.5%
	Doc generation	44	5.57.10 7	18	17	-5.5%
	Arcade Game	61	3.3.10 9	25	14	-44.0%
	Model Transformation 88	1.65.10 13	40	26	-35.0%
	Coche Ecologico	94	2.32.10 7	114	92	-19.3%
	Electronic shopping	287	2.26.10 49	62	37	-40.3%
	Average	74.4	-	34.7	25	-29.2 %
			Table I
	NUMBER OF TEST CONFIGURATIONS ENFORCING PAIRWISE

http://www.irisa.fr/lande/gotlieb/resources/Pacogen/Pacogen.html

http://www.splot-research.org/

www.irisa.fr/lande/gotlieb/resources/Pacogen/


Matthieu Carlier
email: matthieu.carlier@inria.fr

Arnaud Gotlieb

Filtering by ULP maximum

Constraint solving over floating-point numbers is an emerging topic that found interesting applications in software analysis and testing. Even for IEEE-754 compliant programs, correct reasoning over floating-point computations is challenging and requires dedicated constraint solving approaches to be developed. Recent advances indicate that numerical properties of floating-point numbers can be used to efficiently prune the search space. In this paper, we reformulate the Marre and Michel property over floating-point addition/subtraction constraint to ease its implementation in real-world floating-point constraint solvers. We also generalize the property to the case of multiplication/division in order to benefit from its improvements in more cases.

I. INTRODUCTION

Numerical programs used in critical software systems must be thoroughly tested before being embedded within a final product. In the presence of floating-point computations, this involves 1) producing test inputs [START_REF] Miller | Automatic generation of floatingpoint test data[END_REF] targeted to reveal faults and 2) predicting the program results expected on those computations. For this latter problem, called the oracle problem [START_REF] Weyuker | On testing non-testable programs[END_REF], several techniques have been proposed, such as using multiple program executions to check the results in the data diversity approach [START_REF] Ammann | Data diversity: An approach to software fault tolerance[END_REF], or using the Abstract Interpretation framework [START_REF] Cousot | Abstract interpretation : A unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] to estimate the deviance of the floating-point results w.r.t. an interpretation over the reals [START_REF] Goubault | Static analyses of the precision of floating-point operations[END_REF], or else using perturbation techniques to evaluate the stability of a numerical program [START_REF] Tang | Perturbing numerical calculations for statistical analysis of floating-point program (in)stability[END_REF]. Constraint Programming techniques have also been used to tackle the oracle problem: [START_REF] Lozano | Testing continuous double auctions with a constraint-based oracle[END_REF] proposed a Gecode model to evaluate the expected results of a trading system's continuous double auction, [START_REF] Gotlieb | Automated metamorphic testing[END_REF] proposed using the SICStus Prolog clpfd library to generate test inputs that violate properties over imperative programs, etc. However, in the presence of floating-point computations, very few approaches addressed the former problem, namely the test inputs generation problem. Thirty years ago, Miller and Spooner [START_REF] Miller | Automatic generation of floatingpoint test data[END_REF] proposed to automatically find a floatingpoint test input that exercises a given program path, by minimizing a cost function which evaluates the distance between the currently executed path and the expected one. Their work opened the door for search-based test data generation methods [START_REF] Korel | Automated software test data generation[END_REF], [START_REF] Mcminn | Search-based software test data generation: A survey[END_REF] that share great similarities with local search techniques in Constraint Programming [START_REF] Arcuri | Theoretical analysis of local search in software testing[END_REF]. However, these approaches are only based on program executions and do not rely on symbolic reasoning. Thus, they cannot be used to study path feasibility, i.e. to decide whether a possible execution path is feasible or not in the program. In addition, these techniques can be stuck in local minima without being able to provide a meaningful result [START_REF] Arcuri | Theoretical analysis of local search in software testing[END_REF]. An approach to tackle these problems combines both the program execution and the symbolic reasoning [START_REF] Godefroid | Dart: directed automated random testing[END_REF], but this kind of reasoning over floating-point computations is hard and requires dedicated filtering algorithms [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF], [START_REF] Michel | Solving constraints over floating-point numbers[END_REF]. In summary, solving constraints over floating-point numbers allows us to generate test inputs that exercise a selected behaviour of the program under test. This approach is currently implemented in three solvers: the FPCS solver [START_REF] Blanc | The v3f project[END_REF], FPSE for C programs [START_REF] Botella | Symbolic execution of floating-point computations[END_REF] which is available for experiments 1and Gatel, the test data generator for Lustre programs [START_REF] Blanc | The v3f project[END_REF], [START_REF] Marre | Test selection strategies for lustre descriptions in gatel[END_REF].

A promising approach to improve the filtering capabilities of constraints over floating-point variables consists in using numerical properties of floating-point computations. For linear constraints, this led to a relaxation technique where floating-point numbers and constraints are converted into constraints over the reals by using Linear Programming approaches [START_REF] Rueher | Approximating floating-point operations to verify numerical programs[END_REF]. For interval-based consistency approaches, Marre and Michel proposed in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] to exploit the floatingpoint representation in filtering algorithms for both the addition and subtraction constraints.

This paper is concerned with the reformulation of the Marre and Michel property in terms of filtering by maximum ULP (Units in the Last Place), in order to ease its implementation in real-world constraint solvers such as FPSE or FPCS. In addition, we generalize the property to the multiplication and division operators for benefiting of its improvements in more cases. This generalization has an interesting impact on non-linear problems where no other techniques yet exist to improve the search space filtering.

The rest of the paper is organized as follows. Next section briefly presents the IEEE-754 standard of binary floatingpoint numbers and introduces the notations used throughout the paper. Section 3 recalls the basic principles of intervalbased consistency techniques over floating-point variables and constraints. Section 4 presents our generalization of the Marre and Michel property, while section 5 concludes the paper.

II. PRELIMINARIES

A. IEEE-754

This section introduces the arithmetical model specified by the IEEE-754 standard for binary floating-point arithmetic [START_REF]Standard for binary floating-point arithmetic[END_REF].

IEEE-754 specifies two basic binary floating-point formats (single and double) and two extended formats. A floating-point number is noted (-1) s a.m × 2 e where s is the sign bit, a is the hidden bit, m is the significand and e is the biased exponent. The single format occupies 32 bits (1 bit for the sign, 8 for the exponent and 23 for the significand) while the double occupies 64 bits (1 bit for the sign, 11 for the exponent and 52 for the significand). Each format defines several classes of numbers: normalized numbers, denormalized numbers, signed zeros, infinities and NaNs (which stands for Not-a-Number). For the single format, normalized numbers corresponds to an exponent value 0 < e < 255 and a value given by the formula: (-1) s 1.m 2 e-127 . Denormalized numbers correspond to an exponent e = 0 and a value given by (-1) s 0.m 2 -126 where m = 0. Note that the significand possesses a hidden bit which is 1 for normalized numbers and 0 for denormalized. Note also that the bias is equal to 127 for the single format and the exponent is -126 for denormalized numbers. There are two infinities (noted +INF, -INF with e = 255, m = 0) and two signed zeros (noted +0.0, -0.0 with e = 0, m = 0). NaNs (e = 255, m = 0) are used to represent the results of invalid computations such as a division or a subtraction of two infinities. They allow the program execution to continue without being halted by an exception. IEEE-754 indicates five types of rounding directions: toward negative infinity (down), toward positive infinity (up), toward zero (chop) and toward the nearest representable value, with two flavors, tail-to-even or tail-toaway in which respectively values with even mantissa are preferred or values away from zero. The to-the-nearest tailto-even value of a real x will be noted •(x). All rounding modes are monotonic, i.e., ∀x y

∈ R, x ≤ y ⇒ •(x) ≤ •(y).
The most important requirement of IEEE-754 arithmetic is the accuracy of floating-point computations: each of the following operations, add, subtract, multiply, divide, square root, remainder, conversions and comparisons, must deliver to its destination the exact result if possible or the floatingpoint number that requires the least modification of the exact result w.r.t. the prescribed rounding mode and the result format destination. It is said that these operations are correctly rounded For example, the single-format result of 999999995904 + 10000 is2 999999995904 which is the single-format floating-point number nearest to the exact result over the reals. This example shows that the accuracy requirement of IEEE-754 does not prevent surprising results from arising (the second operand is absorbed by the addition operator).

B. Notations

R denotes the set of reals while F denotes an idealized finite set of numerical binary floating-point numbers, defined from a given IEEE-754 format. Throughout the paper, we will consider only floats having a numerical binary representation in the single or the double format, excluding denormalized numbers and NaNs, but including -INF, +INF and zeros. Considering this idealized set is advantageous to simplify the properties and avoid too technical details about denormalized numbers. But, extensions are mentioned when available. A real decimal constant (such as 1.0e12) denotes a floating-point value, and thus, has to be understood as its nearest floating-point number (i.e., as 999999995904). Henceforth x + (resp. x -) denotes the smallest (resp. greatest) floating-point number strictly greater (resp. smaller) than x w.r.t. the considered IEEE-754 format. We denote min the exponent of the smallest normalized numbers in absolute values. Thus, the smallest positive normalized number is 1.0 . . . 0 × 2 min . In our idealized set of numbers, min is the smallest possible exponent for a float. On the opposite, max denotes the greatest possible exponent for a float distinct from -INF or +INF. f max is the greatest representable numerical float. Its value is 1.1 . . . 1 × 2 max and we have f + max = +INF. Arithmetical operations over the floats will be noted using the four operators: ⊕, ⊖, ⊗ and ⊘, corresponding respectively to +, -, * , / over the reals. According to IEEE-754, they are defined with the rounding operator • as follows:

x ⊕ y = •(x + y), x ⊖ y = •(x -y), x ⊗ y = •(x * y), x ⊘ y = •(x/y)
⊙ denotes any of ⊕, ⊖, ⊗ and ⊘. A floating-point variable x will be associated an interval of possible floating-point values, noted x ∈ [x, x] where x denotes the smallest float of x and x its greatest value and x ≤ x. Finally, mid(a, b) denotes the floating-point number at the middle of a and b, which may be a floating-point number of a wider format than of its operands.

III. BACKGROUND ON CONSTRAINT SOLVING OVER FLOATING-POINT VARIABLES

A. Interval-based consistency on arithmetical constraints

[4], [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF] contain formulas for projectors in a intervalbased consistency approach to constraint solving over the floats. It is worth distinguishing direct from indirect projectors as constraints come from program analysis. Roughly speaking, when analyzing an imperative program, assignments are considered as if they were equality constraints, where the assigned variable is labelled with a fresh name. For example, the assignment i++; is translated into the

Addition : z = x ⊕ y z = x ⊕ y, z = x ⊕ y (direct) x = mid(z, z + ) ⊖ y (1 st indirect) x = mid(z, z -) ⊖ y y = mid(z, z + ) ⊖ x (2 nd indirect) y = mid(z, z -) ⊖ x Substraction : z = x ⊖ y z = x ⊖ y, z = x ⊖ y (direct) x = mid(z, z + ) ⊕ y (1 st indirect) x = mid(z, z -) ⊕ y y = x ⊖ mid(z, z -) (2 nd indirect) y = x ⊖ mid(z, z + )
Figure 1. Formulas for direct/indirect projectors equality constraint i 2 = i 1 + 1. Hence, there are two distinct kinds of projection. The projection over variable i 2 is called direct projection while the projector over i 1 is called indirect projection. When solving constraints over floating-point variables, the formulas for direct and indirect projections may be different in order to improve the filtering results. Fig. 1 recalls the formulas used for implementing the interval-based addition/subtraction projectors, while (nonoptimal) formulas for product/division can be found in [START_REF] Botella | Symbolic execution of floating-point computations[END_REF], [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF].

B. The Marre and Michel property

This section presents the Marre and Michel property published in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] for improving the filtering of the addition/subtraction projectors. The idea behind this property comes from the representation of floating-point numbers among the reals: the greater a float is, the greater the distance between it and its immediate successor is. More precisely, for a given float x with exponent n, if ∆ = x + -x, then for y of exponent n + 1 we have y + -y = 2 * ∆.

Fig. 2 gives an intuitive view of the property for subtraction. Let z = y ⊖ x and suppose that z ∈ [v z , v z ] is a strictly positive constant, then the Marre and Michel property says that it exists two greatest values y m and x m such that y m ⊖ x m = v z . It is worth noticing that y m and x m depends neither on y nor on x and do not necessarily belong to them. For example, the existence of x m (similar for y m ) can be explained on Fig. 2

as follows. Consider x

+ m , then for all floating-point v, v ≤ y m =⇒ v -x + m ≤ A and y m < v =⇒ B < v -x + m with A < v z < B. Thus, x +
m cannot be part of the solution. Considering values greater than x + m leads to equivalent inequalities with some A ′ , B ′ such that A ′ ≤ A and B ≤ B ′ . property 3.1: [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] Let a variable z ∈ [z, z] such that 0 < z < z, if z = x ⊖ y, then upper bounds of x and y can be computed using the following formulas. Let ζ be a normalized floating-point number of z such that

ζ =          1.0 . . . 0 × 2 e z iff ez = e z z iff z = 1.0 . . . 0 × 2 ez 1.b 2 . . . b i+1 0 . . . × 2 e z iff z = b 1 .b 2 . . . b i 0b i+2 . . . bp × 2 n z = b 1 .b 2 . . . b i 1b ′ i+2 . . . b ′ p × 2 n with b i+2 = b ′ i+2
and let nb z be the number of zeros in the significand of ζ. There do not exist any value x ′ , strictly greater than β and y ′ , strictly greater than α such that x ′ ⊖ y ′ = ζ where

α = 1.1 . . . 1 × 2 e ζ +nbz β = α ⊕ ζ
The above theorem can be extended to the case where ζ is a denormalized floating-point number by considering that the hidden bit is 0 and by shifting the significand to its first non-zero bit. In the case where z contains only strictly non-negative numbers, the above theorem can be used as is for improving the filtering algorithm of both the addition/subtraction constraints. When z contains only strictly negative numbers, an argument related to the symmetry of floating-point representation can be used as well. For similar reasons, the Marre and Michel property can be used for the filtering algorithm of both lower bounds of variables x and y. However, the above theorem cannot be extended to an interval of z that contains zero. Fig. 3 summarizes the formulas obtained with the Marre and Michel property over floating-point addition/subtraction constraints after exploiting symmetries.

IV. FILTERING BY MAXIMUM ULP

This section reformulates the Marre and Michel property by considering the properties that should be verified by a function δ ⊙ to deduce an optimal interval-consistency based filtering algorithm for ⊙. It also generalizes the property to product/division projectors. In this paper, the filtering algorithms that result from this generalization are collectively called filtering by maximum ULP to refer to the maximal distance between two successive floats.

A. Upper bound

Let δ ⊙ : F → F + be a function that satisfies the following properties:

∀z ∈ F , z = 0 ⇒ ∃y, δ ⊙ (z) ⊙ y = z ∀z z ′ ∈ F , z = 0 ⇒ z ′ > δ ⊙ (z) ⇒ ∄y, z ′ ⊙ y = z
Roughly speaking, those properties say that δ ⊙ returns the smallest float that permits one to recover z. The following is entailed by both properties: property 4.1: Let {v 1 , . . . , v n } be a set of non-zero floats and m a float such that ∀i, δ ⊙ (m) ≥ δ ⊙ (v i ). Then, for all floats m ′ > δ ⊙ (m), there does not exist a float y such that As a consequence, by considering the constraint z = x ⊙ y where 0 / ∈ z, we got that the smallest float m greater than the upper bound of the co-domain of δ ⊙ on z (i.e.,

m ′ ⊙ y = v i . v z x m x + m y m y + m v z ∆ 2 * ∆ B A 0 A = y + m -x + m B = y m -x + m Figure 2. Existence of a limit such that x + -x ≥ vz z > 0 z < 0 z = x ⊖ y x ∈ [-α, β], y ∈ [-β, α] x ∈ [-β, α], y ∈ [-α, β] z = x ⊕ y x ∈ [-α, β], y ∈ [-α, β] x ∈ [-α, β], y ∈ [-α, β]
∀v z ∈ z, δ ⊙ (m) ≥ δ ⊙ (v z )) such that δ ⊙ (m) is the upper bound of x. It is worth says that m belongs to z.
Lastly, if we have a function δ ′ ⊙ which satisfies all the above mentioned properties but where the operands of ⊙ are inverted in the conclusions, i.e.,

∀z ∈ F , z = 0 ⇒ ∃x, x ⊙ δ ′ ⊙ (z) = z ∀z z ′ ∈ F , z = 0 ⇒ z ′ > δ ′ ⊙ (z) ⇒ ∄x, x ⊙ z ′ = z then δ ′
⊙ permits one to estimate the upper bound of y in the above constraint. Moreover, when ⊙ is a commutative operator (e.g., ⊕, ⊗), the properties of δ ⊙ and δ ′ ⊙ are equivalent and δ ⊙ = δ ′ ⊙ .

B. Lower bound

For computing the lower bound, we consider the function δ ⊙ : F → F -having the following properties:

∀z ∈ F , z = 0 ⇒ ∃y, δ ⊙ (z) ⊙ y = z ∀z z ′ ∈ F , z = 0 ⇒ z ′ < δ ⊙ (z) ⇒ ∄y, z ′ ⊙ y = z
As for upper bound, these two properties entailed a property similar to property 4.1. Hence, considering the constraint z = x ⊙ y, the value of z which minimizes δ ⊙ permits to deduce the lower bound of x. example 4.1: Consider the following constraints,

x ∈ [-1.0 × 2 50 , 1.0 × 2 50 ] y ∈ [-1.0 × 2 30 , 1.0 × 2 30 ] z ∈ [1.0, 2.0] z =
x ⊕ y we obtain the following results after filtering without maximum ULP

x = mid(z, z + ) ⊖ y = 1.0 × 2 30 x = mid(z, z -) ⊖ y = -1.0 × 2 30 y = mid(z, z + ) ⊖ x = 1.0 × 2 50 y = mid(z, z -) ⊖ x = -1.0 × 2 50 using maximum ULP ζ = 1.0 . . . 0 × 2 1 and then nb z = 23 α = 1.1 . . . 1 × 2 1+nbz β = α ⊕ ζ = 1.0 . . . 0 × 2 25 x = β = 1.0 . . . 0 × 2 25 x = -α = -1.1 . . . 1 × 2 24 y = β = 1.0 . . . 0 × 2 25 y = -α = -1.1 . . . 1 × 2 24
This example shows that the Marre and Michel property (called filtering by maximum ULP, see below) permits one to get tighter results than classical interval-consistency based filtering. Note that other trivial examples show that the opposite may also arise and therefore it is worth computing the intersection of both filtering to get optimal results. Other examples that show those phenomena can be found in [START_REF] Carlier | Anr u3cat project report -wp1: Floating-point computations[END_REF].

C. Filtering by maximum ULP on addition/subtraction

This section introduces only the functions δ ⊕ and δ ⊕ . The functions δ ⊖ and δ ⊖ can be deduced from symmetries, as explained in Sec. III-B and [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF].

definition 4.1: The function δ ⊕ is defined as follows (i and n are defined such |z| = 1.z 1 . . . z i 10 . . . 0 × 2 n ):

δ ⊕ (0) = 0 δ ⊕ (z) = z + 1.1 . . . 1 × 2 n+(p-i-1) for z > 0 δ ⊕ (z) = 1.1 . . . 1 × 2 n+(p-i-1)
for z < 0 which is always a float as z + 1.1 . . . 1 × 2 n+(p-i-1) = 0.0. . .01z 1 . . .z i 10. . .0 × 2 n+(p-i) ⊕ 0.1. . .11 1 . . . 1 1 × 2 n+(p-i) 1.0. . .01z 1 . . .z i 00. . .0 × 2 n+(p-i) δ ⊕ respects both properties stating for all non-zero z, δ ⊕ (z) returns the smallest float that permits to recover z. δ ⊕ is defined as δ ⊕ (z) = -δ ⊕ (-z) and its properties are entailed by properties of δ ⊕ . definition 4.2: For a variable z of normalized floatingpoint values such that 0 does not belong to z domain, the value in z maximizing δ ⊕ (and minimizing δ ⊕ ) is defined by destr(z) as: Definition 4.2 can easily be extended to intervals with denormalized numbers and intervals composed of negative floating-point values. But it is useless to extend it with intervals containing both negative and positive values as such an interval contains 0 and the property is not applicable. As a result, we got formulas that compute filtering similar to those resulting from the Marre and Michel property. So, the filtering by maximum ULP presented here is similar to the filtering of [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] which is restricted to addition/subtraction operators.

D. Filtering by maximum ULP on the product

Fig. 4 gives an intuitive view of ULP maximum over the floating-point multiplication. Let z = x ⊗ y and suppose that z ∈ [v z , v z ] is a strictly positive constant, then there exists a greatest values x m such that, its exists a value y, x m ⊗ y = v z . As for property 3.1, this value depends neither on x nor on y and do not necessarily belong to them. Let explain the existence of x m on Fig. 2. Consider x m such that x m ⊗ 0 + = v z 3 , then by monotonicity of ⊗, v z < x + m ⊗ 0 + (B in the figure). Thus, still by monotonicity, there is no strictly positive value y such that x + m ⊗ y = v z . By sign analysis, there is no negative or nil value y such that x + m ⊗ y = v z . δ ⊗ (z) is defined as δ ⊗ (z) = |z| * 2 p-min and δ ⊗ (z) = -δ ⊗ (z). The value of an interval z which maximizes δ ⊗ (respecting minimizes δ ⊗ ) is the one that has the greatest absolute value. Hence, for v z an element of z (which does not contain 0) that maximizes δ ⊗ then δ ⊗ (v z ) (resp. δ ⊗ (z)) is an upper bound (resp. a lower bound) of x w.r.t. the constraint z = x ⊗ y. As the product is commutative, the same function can be used for y. Details can be found in [START_REF] Carlier | Anr u3cat project report -wp1: Floating-point computations[END_REF]. As an example, consider the following problem: example 4.2:

x

∈ [-INF, +INF], y ∈ [-INF, +INF], z ∈ [0 + , 1.0 × 2 -30 ], z = x ⊗ y
We got δ ⊗ (1.0 . . . 0 × 2 -30 ) = 1.0 . . . 0 × 2 -30 * 2 23-(-127)-1 = 1.0 . . . 0 × 2 119 and the filtering of the domain of x and y yield:

x ∈ [-1.0 . . . 0 × 2 119 , 1.0 . . . 0 × 2 119 ], y ∈ [-1.0 . . . 0 × 2 119 , 1.0 . . . 0 × 2 119 ]
3 this value always exists since this multiplication is equivalent to an exponent shifting This filtering only applies when the value returned by δ ⊗ is smaller than the greatest representable float distinct from +INF (i.e., when z < 1 × 2 -p ). However, this filtering is useful when zero knowledge is available on the domains of x or y such as in the above example.

E. Filtering by maximum ULP on division

δ ⊘ (z) is defined as δ ⊘ (z) = |z| ⊗ f max . For filtering x in constraint z = x ⊘ y, it suffices to take the value v z that maximize δ ⊘ . This value is max(|z|, |z|). By considering δ ⊘ (z) = -δ ⊘ (z), it is also possible to filter lower bounds.

F. Synthesis

In the previous sections, the following functions have been defined:

V. CONCLUSION

This paper is concerned with constraint solving over floating-point computations. Filtering consistency based algorithms are currently the privileged technology for attacking this problem. This paper reformulates the Marre and Michel property proposed in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] for improving the filtering capabilities of the addition and subtraction operators. It also proposes to generalize the property to the case of multiplication and division, something that has not been proposed elsewhere. Further work includes the exploration of other properties based on linearization of floating-point computations, such as those proposed in [START_REF] Rueher | Approximating floating-point operations to verify numerical programs[END_REF].

0 v z B x m x + m B = x + m ⊗ 0 + v z = x m ⊗ 0 +

Figure 3 .

 3 Figure 3. Summary of the filtering formulas for the Marre and Michel property

  destr(z) = 1.0 . . . 0 × 2 n if e z = e z with n = e z destr(z) = 1.b 2 . . . b i a0 . . . × 2 n if z = 1.b 2 . . . b i b i+1 . . . × 2 n z = 1.b 2 . . . b i b ′ i+1 . . . × 2 n with b i+1 = b ′ i+1 and a = 0 if 1.b 2 . . . b i 0 . . . 0 × 2 n = z 1 otherwise

  δ ⊕ (z) = z + 1.1 . . . 1 × 2 n+(p-i-1) if z = +1.z 1 . . . z i 10 . . . 0 × 2 n δ ⊕ (z) = 1.1 . . . 1 × 2 n+(p-i-1) if z = -1.z 1 . . . z i 10 . . . 0 × 2 n δ ⊕ (z) =δ ⊕ (-z) δ ⊗ (z) = |z| * 2 p-min δ ⊗ (z) =δ ⊗ (z) δ ⊘ (z) = |z| ⊗ fmax (for denormalized z) δ ⊘ (z) =δ ⊘ (z)For each constraint, when 0 ∈ z, interval consistency based filtering yields toConstraint x ⊆ y ⊆ z = x ⊕ y(z > 0) [ δ ⊕ (ζ ), δ ⊕ (ζ )] [ δ ⊕ (ζ ), δ ⊕ (ζ )] z = x ⊕ y(z < 0) [δ ⊕ (ζ ′ ),δ ⊕ (ζ ′ )] [δ ⊕ (ζ ′ ),δ ⊕ (ζ ′ )] z = x ⊖ y(z > 0) [ δ ⊕ (ζ ), δ ⊕ (ζ )] [δ ⊕ (ζ ),δ ⊕ (ζ )] z = x ⊖ y(z < 0) [δ ⊕ (ζ ′ ),δ ⊕ (ζ ′ )] [ δ ⊕ (ζ ′ ), δ ⊕ (ζ ′ )] z = x ⊗ y [δ ⊗ (M ), δ ⊗ (M )] [δ ⊗ (M ), δ ⊗ (M )] z = x ⊘ y [δ ⊘ (M ), δ ⊘ (M )]with ζ = destr(z), ζ ′ = destr(-z) and M = max(|z|, |z|).

Figure 4 .

 4 Figure 4. Existence of a limit such that ∀y, x ⊗ y = vz

http://www.irisa.fr/celtique/carlier/fpse.html

These numbers can be exactly represented by single binary FP numbers.


Mounira Bachir
email: mounira.bachir@inria.fr

Frederic Brault
email: frederic.brault@inria.fr

Sid-Ahmed-Ali Touati
email: sid.touati@uvsq.fr

Albert Cohen
email: albert.cohen@inria.fr

Loop Unrolling Minimisation in the Presence of Multiple Register Types: a Viable Alternative to Modulo Variable Expansion

. However, the degree of unrolling should be minimised to control code size and hence I-cache performance.

. In this paper, we complete our study with an improved algorithm for minimising kernel loop unrolling resulting from cyclic register allocation in the presence of multiple register types showing that considering all register types in conjunction provides a lower unrolling degree than considering each register type in isolation. In addition, we integrate our solution within a real world embedded system compiler: st200cc for the ST2xx family of VLIW embedded processors and compare it to MVE. Our large set of experiments on both high performance and embedded benchmarks (SPEC2000, SPEC2006, MEDIABENCH and FFMPEG) demonstrates the practical applicability and the benefits of our approach.

Introduction

Many high performance or soft real-time applications, such as telecom and image processing, exhibit intensive computations in loops. Software pipelining (SWP) is an important instruction scheduling technique for improving the execution rate of these applications, by exploiting the instruction-level parallelism (ILP) in loops [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF][START_REF] Rau | Code generation schema for modulo scheduled loops[END_REF].

When a loop is software pipelined, we cannot use regular register allocation algorithms because of self-interferences in the usual interference graph [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF][START_REF] Lelait | Contribution à l'allocation de registres dans les boucles[END_REF][START_REF] Fisher | Embedded Computing: a VLIW Approach to Architecture, Compilers and Tools[END_REF]. In compiler con-struction, when no hardware support is available, kernel loop unrolling is the method of code generation that does not alter the initiation interval after software pipelining. In fact, unrolling the loop allows us to avoid introducing unnecessary move and spill operations after a periodic register allocation.

Our objective in this research effort is somewhat different: we are interested in the minimal loop unrolling factor which allows a periodic register allocation for software pipelined loops in the presence of multiple register types. We also focus on generating an embedded VLIW loop with reduced spill and code size, without necessarily reducing the execution speed. Unlike high performance and intensive computational fields, we are asked to generate a code with equivalent performance but with less memory operations. We believe that an important code quality criterion is to have a reduced amount of memory requests upon the condition of not altering ILP scheduling and execution speed if possible.

In a target architecture with multiple register types (a.k.a. classes), some of state-of-the-art algorithms [START_REF] Touati | Early periodic register allocation on ilp processors[END_REF][START_REF] Lelait | Contribution à l'allocation de registres dans les boucles[END_REF] propose to compute the sufficient unrolling degree that we should apply to the loop so that it is always possible to allocate the variables of each register type with a minimal number of registers (MAXLIVE [START_REF] Huff | Lifetime-sensitive modulo scheduling[END_REF]). However, periodic register allocation has long suffered from an important drawback: the resulting unrolling factor can be very high. We recently proposed a solution to dramatically reduce this cost when considering a single register type [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF].

In this article, we extend our previous results in two directions:

1. We first provide an efficient algorithm for minimising loop unrolling in the presence of multiple register types, for instance, T = {int, f loat, branch}. The minimal loop unrolling degree is reached by exploiting the unused registers to adjust the different weights of reuse circuits generated for each register type, looking for a good distribution of those registers concurrently over all types.

Note that, the mathematical nature of the two problems are not equivalent (single register type vs. multiple register types). Indeed, if we have two or three types of registers, then we have two or three optimisation problems to solve [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF]. This paper solves them simultaneously by enumerating the feasible unrolling degrees and searching for the minimal kernel loop unrolling among them.

2. Second, we implemented and integrated our code optimisation method within the industrial compiler of STMicroelectronics, allowing us to make extensive experiments on real world high performance and embedded applications.

This paper is organised as follows. Sect. 2 presents some relevant related work on code generation for periodic register allocation. Sect. 3 formalises the problem of minimising the loop unrolling degree in the presence of multiple register types. Sect. 4 details our algorithmic solution for minimising loop unrolling. Sect. 5 presents detailed experimental results on standard benchmarks (MEDIABENCH, SPEC2000, SPEC2006, STMicroelectronics benchmark), showing that our method is efficient in practice. Finally, we conclude.

Related Work

We detail in the following section two techniques that we often refer in this paper. For a comprehensive bibliography on register allocation problems we refer the reader to [START_REF] Touati | Register Pressure in Instruction Level Parallelism[END_REF][START_REF] Lelait | Contribution à l'allocation de registres dans les boucles[END_REF] which contain an extensive list of contributions in the area.

Modulo Variable Expansion

This method, proposed by Lam [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF][START_REF] Rau | Register allocation for software pipelined loops[END_REF], defines a minimal unrolling degree to permit code generation after a given periodic register allocation. This unrolling degree is obtained by dividing the length of the longest lifetime (LT ) of all register types (max t∈T (max v LT v,t )) by the initiation interval maxt∈T (maxv LTv,t) II

. However, this method does not guarantee for each type t, a register allocation with a minimal number of registers equal to MAXLIVE, the number of values simultaneously alive for this register type. In general, as we will see in Sect. 5, it may lead to unnecessary spills or move operations negating the benefits of SWP. These extra spill or move operations may increase the initiation interval of the SWP.

SIRA Reuse Graphs

Reuse graphs are a generalisation of the work previously published by [START_REF] De Werra | On a graph-theoretical model for cyclic register allocation[END_REF][START_REF] Hendren | A register allocation framework based on hierarchical cyclic interval graphs[END_REF] and are used inside a framework called SIRA [START_REF] Touati | Early periodic register allocation on ilp processors[END_REF]. Unlike the previous approaches for periodic register allocation, reuse graphs can be used before or after software pipelining to generate a move-free or a spill-free periodic register allocation. Reuse graphs give a formal method to generate code without spill and without move operations if loop unrolling is applied. A brief overview of reuse graphs is illustrated in Fig. 1. Fig. 1 ) share the same destination register. Hence, reuse graphs completely define a periodic register allocation for a given loop, either before SWP (unscheduled loop) or after SWP (already scheduled loop). A formal theorem in [START_REF] Touati | Early periodic register allocation on ilp processors[END_REF] proves that the number of allocated registers of type t is equal to R min,t = (u,v) ν t u,v if we unroll the loop with a factor equal to α t . Fig. 1(b) represents a reuse graph that allocates 3 + 2 = 5 registers of type t 1 and 3 + 1 + 3 = 7 registers of type t 2 .

Each register type t requires an unrolling factor α t . If the reuse graph G t r contains multiple reuse circuits C 1 , • • • , C k , then the weight of each reuse circuit is defined by µ i,t = (u t ,v t )∈Ci ν t u,v . The unrolling degree of type t is then equal to α t = lcm(µ 1,t , • • • , µ k,t ). For instance, the unrolling degree of G t2 r in Fig. 1(b) is equal to α t2 = lcm(3 + 1, 3) = lcm(4, 3) = 12. Similarly, α t1 = lcm(3, 2) = 6. The global unrolling degree that is valid for all register types concurrently is equal to α = lcm t∈T (α t ). For Fig. 1(b), α = lcm(α t1 , α t2 ) = lcm(6, 12) = 12.

(1,0) [START_REF] Lelait | Contribution à l'allocation de registres dans les boucles[END_REF][START_REF] Touati | Early periodic register allocation on ilp processors[END_REF] (3, 1)

(1,1) (2,0) (1,2) (a) Initial DDG (b)
Reuse Graphs for Register Types t1 and t2

u 3 u 4 u 2 u 1 u 2 u 4 Register type t 2 u 3 u 1 ν u1,u1 = 3 ν u4,u2 = 1 ν u2,u4 = 3 ν u1,u1 = 3 ν u3,u3 = 2 u 1
Register type t 1

Loop Unrolling Problem

Code generation methods using loop unrolling [START_REF] Touati | Early periodic register allocation on ilp processors[END_REF][START_REF] Lelait | Contribution à l'allocation de registres dans les boucles[END_REF] have long suffered from a drawback: the resulting unrolling factor α can be high. We recently proposed a solution to minimise the unrolling degree when considering a single register type [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF]. However, in the presence of multiple register types, minimising the loop unrolling degree of each type seperately does not lead to the minimal loop unrolling degree for the whole loop. Fig. 2 illustrates an example. We want to minimise the loop unrolling degree of the initial reuse graph in Fig 1(b), where two register types t 1 , t 2 are considered. The initial kernel loop unrolling degree α = 12 is the LCM of α t1 = 6 and α t2 = 12 which are respectively the LCM of the different reuse circuits weights for each register type. In this configuration, let assume that we have R hw,tj = 8 architectural registers in the processor for each register type t j . Hence we have R t1 = R hw,t1 -R min,t1 = 8 -5 = 3 (resp R t2 = 1) remaining registers for register type t 1 (resp t 2 ). By applying the loop unrolling minimisation for each register type separately [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF] (Fig 2(a)), the minimal loop unrolling degree for each register type becomes: α * t1 = 3 for register type t 1 and α * t2 = 4 for register type t 2 . However, the global kernel loop unrolling degree is not minimal α ′ = lcm(α * t1 , α * t2 ) = 12. This paper describes how to find the minimal loop unrolling degree α * for all register types concurrently. Our goal is to exploit the remaining registers of each register type, looking for a good distribution of those registers over all types. In Fig 2(b), the final loop unrolling degree found with this new method called Loop Unrolling Minimisation (LUM) is α * = 4 < α ′ . The minimal number of registers added to each reuse circuit of each type are: r 1,t1 = 1, r 2,t1 = 0, r 1,t2 = 1, r 2,t2 = 0. Note that r i,tj is the number of registers added to the i th reuse circuit of the type t j . LUM guarantees that the new number of allocated registers will not exceed the number of architectural registers for each register type t j ; R alloc,tj ≤ R hw,tj . Minimal Kernel Loop Unrolling:

C2,t 1 = {u3} µ2,t 1 = 2 + 1 µ1,t 1 = 3 C1,t 2 = {u1} µ1,t 2 = 3 + 1 µ2,t 2 = 4 u 2 u 4 u 1 u 3 u 1 C2,t 2 = {u2, u4} Register type t1 : Rhw,t 1 = 8 α * t1 = lcm(3, 3) = 3 Ralloc,t 2 = 4 + 4 ≤ Rhw,t 2 α * t2 = lcm(4, 4) = 4 Ralloc,t 1 = 3 + 3 ≤ Rhw,t 1 C1,t 1 = {u1} C2,t 1 = {u3} µ2,t 1 = 2 C1,t 2 = {u1} µ1,t 2 = 3 + 1 u 2 u 1 u 3 u 1 µ2,t 2 = 4 µ1,t 1 = 3 + 1 u 4 α * = lcm(α ′ t1, α ′ t2) = 4 C1,t 1 = {u1}
Register type t2 : Rhw,t 2 = 8

Register type t1 : Rhw,t 1 = 8 Register type t2 : Rhw,t 2 = 8

α ′ = lcm(α * t1, α * t2) = 12 α ′ t2 = lcm(4, 4) = 4 Ralloc,t 2 = 4 + 4 ≤ Rhw,t 2 α ′ t1 = lcm(4, 2) = 4 C2,t 2 = {u2, u4} Ralloc,t 1 = 4 + 2 = 6 ≤ Rhw,t 1

Figure 2. Modifying Reuse Graphs to Minimise Loop Unrolling Factor

In the next section, we give a formal description for minimising loop unrolling in the presence of multiple register types.

Loop Unrolling Minimisation Problem

PROBLEM 1 (LUM). Let α be the initial loop unrolling degree and let T = {t 1 , . . . , t n } be the set of register types. For each register type t j ∈ T , let R tj ∈ N be the number of remaining registers after a periodic register allocation for this register type. Let k j be the number of generated reuse circuits. We assume µ i,tj ∈ N is the weight of the i th reuse circuit of the register type t j . For each reuse circuit i and each register type t j , we must compute the additional registers r i,tj such that we find a new periodic register allocation with a minimal loop unrolling degree. This can be described by the following constraints:

1. α * = lcm(lcm(µ 1,t1 + r 1,t1 , . . . , µ k1,t1 + r k1,t1 ), . . . , lcm(µ 1,tn + r 1,tn , . . . , µ kn,tn + r kn,tn )) is minimal

2. ∀t j ∈ T , kj i=1 r i,tj ≤ R tj
Problem 1 proposes to exploit the remaining registers of each type to adjust the weights of the different reuse circuits, looking for an optimal distribution of those registers concurrently over all types in order to find a new periodic register allocation with a minimal loop unrolling degree.

The following section defines the search space S for the minimal kernel loop unrolling α * .

Search Space for Minimal Kernel Loop Unrolling

According to LCM properties and to the formulation of Problem 1, the search space S for the minimal kernel loop unrolling α * is bounded. In fact, three cases arise:

Case 1: No remaining registers for all register types

In this case, the initial loop unrolling degree cannot be minimised α * = α.

Case 2: No remaining registers for some register types

Assume that α j is the loop unrolling degree for the register type t j ∈ T . By the way, α = lcm(α 1 , . . . , α n ). We define the subset T ′ which contains all the register types such that they have no remaining registers after periodic register allocation (T ′ ⊂ T such that

T ′ = {t ∈ T | R t = 0}).
In fact, if there are no registers left for these register types, we cannot minimise their loop unrolling degrees [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF]. Therefore, the minimal global loop unrolling degree α * ≥ α j ∀ t j ∈ T ′ . By assuming α ′ = lcm t∈T ′ (α t ), the search space S is defined as follows:

S = {β ∈ N | β is multiple of α ′ ∧ α ′ ≤ β ≤ α}
Here, each value β can be a potential final loop unrolling degree.

Case 3: All register types have some remaining registers

The final loop unrolling factor α * is a multiple of each updated reuse circuit weight (µ i,tj + r i,tj ) with the number of additional registers (r i,tj ) varied from 0 (no added register for this circuit) to R tj (all the remaining registers are added to this circuit).

Furthermore, if we assume that µ kn,tn is the maximum weight of all the different circuits for all register types (µ kn,tn = max tj (max i µ i,tj )) then α * is a multiple of this specific updated circuit (α * is a multiple of (µ kn,tn +r kn,tn ) with 0 ≤ r kn,tn ≤ R tn ). We notice here that any reuse circuit satisfies this later property, but it is preferable to consider the reuse circuit with a maximal weight because it decreases the cardinality of the search space S. Finally the search space S can be stated as follows:

S = {β ∈ N | β is multiple of (µ kn,tn + r kn,tn ) , ∀r kn,tn = 0, R tn ∧ µ kn,tn ≤ β ≤ α}
After describing the set S of all possible values of α * (case 2 and case 3), the minimal kernel loop unrolling α * is defined as follows:

α * = min{β ∈ S|∀t j ∈ T, ∃(r 1,tj , . . . , r kj ,tj ) ∈ N kj such that: β = lcm(lcm(µ 1,t1 +r 1,t1 , . . . , µ k1,t1 +r k1,t1 ), . . . , lcm(µ 1,tj + r 1,tj , . . . , µ kj ,tj +r kj ,tj ), . . . , lcm(µ 1,tn +r 1,tn , . . . , µ kn,tn +

r kn,tn )) ∧ kj i=1 r i,tj ≤ R tj }
Here arises another problem: how to decide if the value β can be a potential new loop unrolling. A proposition for solving this problem is explained in the next section.

Fixed Loop Unrolling Problem

PROBLEM 2 (Fixed Loop Unrolling). Let β ∈ S be a fixed loop unrolling degree and let T = {t 1 , . . . , t n } be the set of register types. β can be a potential new loop unrolling iff we find for each register type t j ∈ T , a minimal distribution of the remaining registers R tj between its reuse circuits (µ i,tj ) such that this new loop unrolling degree β satisfies the following constraints:

1. β = lcm(lcm(µ 1,t1 + r 1,t1 , . . . , µ k1,t1 + r k1,t1 ), . . . , lcm(µ 1,tn + r 1,tn , . . . , µ kn,tn + r kn,tn ))

2. ∀t j ∈ T kj i=1 r i,tj ≤ R tj
In order to determine if β can be the new kernel loop unrolling, we propose to generalise the LCM-Problem solution described in our previous work [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF] for all register types. The constraints in Problem 2 are LCM-Problem constraints which must be checked for all the register types.

In general, the LCM-Problem proposes to add to each reuse circuit µ i,tj of each register type t j , a minimal number of registers r i,tj from the remaining R tj registers such that µ i,tj + r i,tj is the smallest divisor of the fixed loop unrolling β greater or equal to µ i,tj . In this way, if the additional registers, for each register type, do not exceed the number of remaining registers kj i=1 r i,tj ≤ R tj , then β can be the new loop unrolling degree.

However, in the presence of multiple register types, the meaning is slightly different. β is, in fact, the least common multiple of the loop unrolling for all the register types. On the contrary, if we consider each register type separately, β is not necessarily the least commom multiple of its different updated reuse circuits weights, but a multiple of their least common multiple.

The solution of the Fixed Loop Unrolling Problem constitutes the basis of the solution for Loop Unrolling Minimisation Problem explained in the next section.

Solution for Minimal Loop Unrolling

In order to compute the minimal kernel loop unrolling α * , our solution consists of checking if each value β in the search space S can be a solution for the Fixed Loop Unrolling Problem: it is guaranteed that the minimum of all these values is the minimal loop unrolling degree.

Instead of computing all values β of S which satisfy the Fixed Loop Unrolling Problem and finally taking the minimal one, we describe in Fig. 3 an efficient way to find the minimal α * depending on the construction of the set S. Fig. 3 also illustrates the different cases of the construction of the research space S. The value of each node represents a potential new loop unrolling degree and an edge between two nodes a, b (a → b) means that a < b and a dashed edge between two nodes means the order is unknown. The structure of the search space depends on the availablity of the different types of registers :

• Case 1 (no registers left for all register types): no loop unroll minimisation is possible, α * = α.

• Case 2 (no registers left for some register types): α * is multiple of α ′ , we apply the algorithm of Fixed Loop Unrolling Problem to each node of Fig. 3 until we find a solution or until we reach the last node α

• Case 3: some registers left for all register types: we traverse the set S in the same way as described in [START_REF] Bachir | Post-pass periodic register allocation to minimise loop unrolling degree[END_REF].

If we assume that µ = µ kn,tn (maximum weight of all the different circuits for all register types) and R = R tn (remaining registers for the register type t n ) then we traverse the set S by proceeding line by line. In each line, we apply the solution of Fixed Loop Unrolling Problem to each node in turn until we find a value which satisfies the constraints of Problem 2 or until we arrive at the last line where β = α. If the value β of the node i of the line j is a solution for the Fixed Loop Unrolling Problem, then we have two cases: a) If the value of this node is less than the value of the first node of the next line then we are sure that this value is minimal (α * = β). This is because all the remaining nodes are greater than β (by construction of the set S).

b) Otherwise we have found a new value of unrolling degree which is less than the original α. We note this new value α" and we try once again to min-imise it until we find the minimal (case a). The search space becomes smaller (S ′ = {β ∈ N|∀r = 0..R :

β is multiple of (µ + r) ∧ (j + 1) × µ ≤ β ≤ α"})
α β = α ′′ α α ′ 2 * α ′ α * α µ µ + 1 µ + 2 µ + R 2 * µ 2 * (µ + 1) 2 * (µ + 2) 2 * (µ + R) 3 * µ 3 * (µ + 1) 3 * (µ + 2) 3 * (µ + R) α * = α β = α ′′

Experimental Results

In our study, we take the ST231 embedded VLIW processor as a target. It is currently the latest processor of the ST2xx family from STMicroelectroncis. ST231 is an integer 32 bits VLIW processor [START_REF] Faraboschi | Lx: a technology platform for customizable vliw embedded processing[END_REF]. The size of the L1 I-cache and Dcache is 32 KB each. For our experiments, we consider the two register types of ST2xx, T = {general, branch}. The number of architectural registers of each register type is configured as follows: R arch,general = 32 × 32 bits , R arch,branch = 4× 1-bit. Our loop unrolling minimisation method is independent of the technique used for periodic register allocation. Consequently, it can be performed after any periodic register allocation technique. For the purpose of our study, we chose to develop a new version of SIRA [START_REF] Touati | Early periodic register allocation on ilp processors[END_REF] implementing a periodic register optimisation with loop unrolling minimisation in the presence of multiple register types. Then, we integrated our SIRA optimiser inside the st200cc compiler from STMicroelectroncis.

To study the efficiency of our loop unrolling minimisation, we conducted an extensive set of experiments on both high performance and embedded benchmarks: SPEC2000, SPEC2006, MEDIABENCH and ST FFMPEG internal benchmark; FFMPEG is a representative video multimedia application which is used with the ST231. The experiments were performed on close to 9000 software pipelined loops.

First, our experiments show that the run-time of our register allocation followed by loop unrolling minimisation is less than 1 second per loop on average (on a 3.4 GHz Pentium D). So, it is fast enough to be included inside an industrial cross compiler such as st200cc.

The following sections present our experiments. We focus on the benefit of our new method comapred to MVE (in terms of code size, spill reduction, move reduction, and II increase). In order to have an idea of the different results, we chose to graph them as boxplots, which are a convenient way of graphically depicting groups of numerical data through their five-number summaries: the smallest observations, lower quartile (Q1 = 25%), median (Q2 = 50%), upper quartile (Q3 = 75%), and largest observations.

Statistics on Minimal Loop Unrolling Factors

Fig. 4 shows numerous boxplots representing the initial loop unrolling degree and the final loop unrolling degree of the different loops per benchmark application. In each benchmark family (FFMPEG, MEDIABENCH, SPEC2000, SPE2006), we note that the loop unrolling degree is reduced significantly from its initial value to its final value.

To highlight the improvements of our loop unrolling minimisation method, we show in Fig. 5 a boxplot for each benchmark family (FFMPEG, SPEC2000, SPEC2006, ME-DIABENCH). We remark that the final loop unrolling of half of the applications is under 3 and that the final loop unrolling of 75% of applications is less than or equal to 5. This compares favourably with the loop unrolling degrees calculated by minimising each register type in isolation. Here, the final loop unrolling degree of half of the applications is under 5 and the final loop unrolling of 75% of the applications is under 7, the final loop unrolling for the remaining loops can reach 50. These numbers demonstrate the advantage of minimising all register types concurrently.

In order to compare our method with modulo variable expansion (MVE [START_REF] Lam | Software pipelining: an effective scheduling technique for vliw machines[END_REF]), we also plot in Fig. 5 the loop unrolling degree obtained thanks to MVE. The MVE heuristic always finds a smaller unrolling degree. However, MVE does not guarantee a register allocation with a minimal number of registers and in general it may lead to unnecessary spills breaking the benefits of SWP. In other words, MVE may require spill code insertion even if MAXLIVE does not exceed the number of architectural registers. This problem occurs in 254 loops of FFMPEG, 405 loops of MEDIABENCH, 1270 loops of SPEC2000 and 571 loops of SPEC2006.

With our loop unrolling register allocation method, we formally guarantee that we do not need, for each register type, more registers than the number of architectural registers. In addition, a suitable quality criterion is to check if the unrolled loops fit in the I-cache. Fortunately, our experimental results show that using loop unrolling technique, all the unrolled loops also fit in the I-cache, even if our unrolling factors are higher than those computed by MVE.

The next section demonstrates that loop unrolling minimisation is a better choice than MVE due to the reduction in spill code and move operations.

SPEC2006 Final Loop Unrolling

Statistics on Spill Code Decrease and Move Operation Decrease

The spill code decrease is computed over all SWP loops. It is evaluated as the difference between the amount of spill code generated without using SIRA, i.e. with MVE, and the amount of spill code generated when using SIRA, all divided by the total amount of spill code generated without using SIRA ( P Spillnosira-P Spillsira P Spillnosira

). The results show that including the new version of SIRA in the st200cc compiler greatly reduces the amount of spill. In fact, the final counts of spills in SWP loops are: FFMPEG=86, MEDIA-BENCH=100, SPEC2000=240, SPEC2006=111.

In addition, we did statistics on the amount of move operations reduced thanks to our method (on all SWP loops). The decrease in move operations is calculated as the difference between the number of the move operations generated without using SIRA, i.e. with MVE, and the number of move operations generated with the use of SIRA, all divided by the number of move operation generated without using SIRA: erations by 10 -20% : we reduce the amount of move operations by 18.92% for FFMPEG, 9.61% for MEDIABENCH, 14.17% for SPEC2000 and 18.21% for SPEC2006.

The way in which SIRA adds arcs to the DDG before software pipelining may in theory increase the critical circuits (minII) and modify ILP scheduling. This is studied in the next section.

Statistics on II increase

Our cyclic register allocations and loop unrolling minimisation are performed before SWP. It may be argued that introducing arcs inside DDGs before software pipelining would alter the ILP scheduling quality, since extra constraints are added. In practice, this is not the case because the usual software pipelining heuristics are not optimal. This means that even if we introduce additional arcs (either with or without increasing minII), this does not necessarily increase the II. The mean of II increase is resp. 1.56% for FFMPEG, 0.05% for MEDIABENCH, 1.66% for SPEC2000 and 0.09% for SPEC2006. As can be seen, the average of II increase is negligible in all benchmarks.

Conclusion

We presented, in this article, an efficient algorithm to minimise kernel loop unrolling resulting from the periodic register allocation in the presence of multiple register types. Our experiments on the embedded VLIW processor ST231 cover a wide range of high-performance and embedded benchmarks (FFMPEG, MEDIABENCH, SPEC2000, SPEC2006). These experiments demonstrate the practical applicability and the benefits of our approach: as our periodic register allocation is applied before software pipelining, it avoids the uncontrolled generation of spill code. Regarding the initiation intervals, the increase is negligible on average, and always less than 1.6% on average. Furthermore, thanks to loop unrolling minimisation, all the loops fit into the I-cache. The total amount of register move operations is reduced, yet still over 18% for FFMPEG and SPEC2006, 14% for MEDIA-BENCH and 9% for SPEC2000.

If code size is more critical, fitting the loops inside the I-cache may not always be satisfactory as a quality criteria. In this case, we have still some opportunities to reduce loop unrolling factors by exploiting move operations. This is our current research direction that may combine loop unrolling with move insertions.

Figure 1 .

 1 Figure 1. SIRA Reuse Graphs

2 (

 2 b) Minimising Loop Unrolling for all Register Types ConjointlyNew Kernel Loop Unrolling:

Figure 3 .

 3 Figure 3. Loop Unrolling Values in the Search Space S

Figure 4 .

 4 Figure 4. Statistics on Loop Unrolling Minimisation

.Figure 5 .

 5 Figure 5. Loop Unrolling Minimisation versus MVE

  (a) shows an initial DDG with two register types t 1 and t 2 . Statements which involve registers of type t 1 are in dashed circles, and those of type t 2 are in bold circles. Statement u 1 writes two results of distinct types. Flow dependence through registers of type t 1 are in dashed arcs, and those of type t 2 are in bold arcs. We associate a reuse graph G t r with each register type t, see Fig.1(b). A reuse graph G t r of register type t contains only the nodes writing inside registers of type t. These nodes are connected by reuse arcs. The existence of a reuse arc (u t , v t ) of distance ν t u,v means that the two operations u t (i) and v t (i + ν t u,v