Roberto Bagnara
email: roberto.bagnara@bugseng.com

Matthieu Carlier

Roberta Gori
email: gori@di.unipi.it

Arnaud Gotlieb
email: arnaud@simula.no

Bugseng

Symbolic Path-Oriented Test Data Generation for Floating-Point Programs

Verifying critical numerical software involves the generation of test data for floating-point intensive programs. As the symbolic execution of floating-point computations presents significant difficulties, existing approaches usually resort to random or search-based test data generation. However, without symbolic reasoning, it is almost impossible to generate test inputs that execute many paths with floating-point computations. Moreover, constraint solvers over the reals or the rationals do not handle the rounding errors. In this paper, we present a new version of FPSE, a symbolic evaluator for C program paths, that specifically addresses this problem. The tool solves path conditions containing floating-point computations by using correct and precise projection functions. This version of the tool exploits an essential filtering property based on the representation of floating-point numbers that makes it suitable to generate path-oriented test inputs for complex paths characterized by floating-point intensive computations. The paper reviews the key implementation choices in FPSE and the labeling search heuristics we selected to maximize the benefits of enhanced filtering. Our experimental results show that FPSE can generate correct test inputs for selected paths containing several hundreds of iterations and thousands of executable floating-point statements on a standard machine: this is currently outside the scope of any other symbolicexecution test data generator tool.

I. INTRODUCTION

During the last decade, the use of floating-point computations in the design of critical systems has become increasingly acceptable. Even in the civil and military avionics domain, which are among the most critical domains for software, floating-point numbers are now seen as a sufficiently-safe, faster and cheaper alternative to fullycontrolled, implementation-based fixed-point arithmetic.

Acceptance of floating-point computations in the design of critical systems took a long time. In fact, rounding errors are difficult to predict and control, and can lead to catastrophic failures. For instance, during the first Persian Gulf War, the failure of a Patriot missile battery in Dhahran was traced to an accumulating rounding error in the continuous execution of tracking and guidance software, causing the death of several civilians [START_REF] Skeel | Roundoff error and the Patriot missile[END_REF]. A careful analysis of this failure revealed that, even though the rounding error obtained at each step of the floating-point computation was very small, the propagation during a long loop-iterating path could lead to dramatic imprecision.

Adoption of floating-point computations in critical systems involves the use of thorough unit testing procedures that are able to exercise complex chains of floating-point operations. In particular, a popular practice among software engineers in charge of the testing of floating-point-intensive computations consists in executing carefully chosen loopiterating paths in programs. They usually pay more attention to the paths that are most likely to expose the system to unstable numerical computations. For critical systems, a complementary requirement is to demonstrate the infeasibility of certain paths, in order to convince a third-party certification authority that certain unsafe behaviors of the systems cannot be reached. As a consequence, software engineers face two difficult problems:

1) How to accurately predict the expected output of a given floating-point computation?1 2) How to find a test input that is able to exercise a given path, the execution of which depends on the results of floating-point computations, or to guarantee that such a path is infeasible? The first problem has been well addressed in the literature [START_REF] Kuliamin | Standardization and testing of mathematical functions[END_REF] through several techniques, either based on multiple related program executions [START_REF] Ammann | Data diversity: An approach to software fault tolerance[END_REF], [START_REF] Chan | Application of metamorphic testing in numerical analysis[END_REF], or on statically-extracted properties of programs [START_REF] Goubault | Static analyses of the precision of floating-point operations[END_REF], or on 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]. In contrast, the second problem received only little attention. Beyond the seminal work of W. Miller and D. L. Spooner [START_REF] Miller | Automatic generation of floating-point test data[END_REF], who proposed to guide the search of floating-point inputs to execute a selected path, few approaches try to exactly reason about floating-point computations. The work in [START_REF] Miller | Automatic generation of floating-point test data[END_REF] paved the way to the development of search-based test data generation techniques, which consist in searching test inputs by minimizing a cost function, evaluating the distance between the currently executed path and a targeted selected path [START_REF] Korel | Automated software test data generation[END_REF]- [START_REF] Arcuri | Theoretical analysis of local search in software testing[END_REF]. Although these techniques enable quick and efficient coverage of testing criteria such as "all decisions," they are unfortunately sensible to the rounding-errors incurred in the computation of the branch distance [START_REF] Arcuri | Theoretical analysis of local search in software testing[END_REF]. Moreover, searchbased test data generation cannot be used to study path feasibility, i.e., to decide whether a possible execution path involving floating-point computations 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 program execution and symbolic reasoning [START_REF] Godefroid | DART: Directed automated random testing[END_REF]. This kind of reasoning requires solving constraints over floating-point numbers in order to generate test inputs that exercise a selected behavior of the program under test. However, solving floating-point constraints is hard and requires dedicated filtering algorithms [START_REF] Michel | Solving constraints over floating-point numbers[END_REF], [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF]. According to our knowledge, this approach is currently implemented in four solvers only: ECLAIR 2 , FPCS [START_REF] Blanc | The V3F project[END_REF], FPSE 3 [START_REF] Botella | Symbolic execution of floating-point computations[END_REF], and Gatel, a 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 some peculiar numerical properties of floating-point numbers. 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] Belaid | Boosting local consistency algorithms over floating-point numbers[END_REF]. For interval-based consistency approaches, B. Marre and C. Michel identified a property of the representation of floating-point numbers and proposed to exploit it in filtering algorithms for addition and subtraction constraints [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF]. In [START_REF] Carlier | Filtering by ULP maximum[END_REF], a reformulation of the Marre-Michel property in terms of filtering by maximum ULP (Units in the Last Place) was proposed in order to ease its implementation in constraint solvers such as FPSE. In addition, the authors sketched a generalization of the property to multiplication and division constraints. This paper is concerned with this challenge. More precisely, it addresses the question of whether the Marre-Michel property can be useful for the automatic solution of realistic test input generation problems. The contributions of the paper are:

1) The filtering algorithm proposed in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] for addition and subtraction is reformulated and corrected.

2) The plan anticipated in [START_REF] Carlier | Filtering by ULP maximum[END_REF] is brought to completion: a uniform framework is thoroughly defined that generalizes the property identified by Marre and Michel to the case of multiplication and division. 3) Our implementation of filtering by maximum ULP in FPSE is presented and critical design choices (e.g., to avoid slow convergence phenomena) are explained. 4) Experimental results are presented on constraint systems that have been extracted from programs engaging into intensive floating-point computations. These results show that the Marre-Michel property and its 2 http://bugseng.com/products/eclair 3 http://www.irisa.fr/celtique/carlier/fpse.html generalization defined in this paper speed up the test inputs generation process. The rest of the paper is organized as follows. Next section presents the IEEE 754 standard of binary floating-point numbers and introduces the notations used throughout the paper. Section III recalls the basic principles of intervalbased consistency techniques over floating-point variables and constraints. Section IV presents our generalization of the Marre-Michel property, while Section V details our implementation of the property in FPSE. Section VI presents our experimental results and analysis. Section VII discusses related work and Section VIII concludes the paper.

II. PRELIMINARIES

A. IEEE 754

This section recalls the arithmetic model specified by the IEEE 754 standard for binary floating-point arithmetic [START_REF]IEEE Standard for Floating-Point Arithmetic, IEEE Std[END_REF].

IEEE 754 binary floating-point formats are uniquely identified by: p ∈ N, the number of significant digits (precision); e max ∈ N, the maximum exponent; e min ∈ N, the minimum exponent (usually 1 -e max ). The single precision format has p = 24 and e max = 127, the double precision format has p = 53 and e max = 1023 (IEEE 754 also defines extended precision formats). An IEEE 754 floating-point number z has the form (-1) s a.m × 2 e where s is the sign bit, a is the hidden bit, m is the significand and the exponent e is also denoted by e z . Each format defines several classes of numbers: normal numbers, subnormal numbers, signed zeroes, infinities and NaNs (Not a Number). The smallest positive normal floating-point number is f nor min = 1.0 • • • 0 × 2 emin = 2 emin and the largest is 1-p ); normal numbers have the hidden bit a = 1. The non-zero floating-point numbers whose absolute value is less than 2 emin are called subnormals: they always have fewer than p significant digits as their hidden bit is a = 0. Every finite floating-point number is an integral multiple of the smallest subnormal magnitude f min = 0.0 • • • 01 × 2 emin = 2 emin+1-p . There are two infinities, denoted by +∞ and -∞, and two signed zeros, denoted by +0 and -0: they allow some algebraic properties to be maintained [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]. NaNs 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.

f max = 1.1 • • • 1×2 emax = 2 emax (2-2
IEEE 754 defines five rounding directions: toward negative infinity (down), toward positive infinity (up), toward zero (chop) and toward the nearest representable value (near); the latter comes into two flavors: tail-to-even or tail-to-away in which values with even mantissa or values away from zero are preferred, respectively. This paper is only concerned with round-to-nearest, tail-to-even, which is, by far, the most widely used. The round-to-nearest, tail-toeven value of a real number x will be denoted by •(x).

All rounding modes are monotonic; in particular, for each x, y ∈ R, x ≤ y implies •(x) ≤ •(y).

The most important requirement of IEEE 754 arithmetic is the accuracy of floating-point computations: add, subtract, multiply, divide, square root, remainder, conversion and comparison operations must deliver to their destination the exact result rounded as per the rounding mode in effect and the format of the destination. It is said that these operations are "exactly rounded."

The accuracy requirement of IEEE 754 can still surprise the average programmer: for example the single precision, round-to-nearest addition of 999999995904 and 10000 (both numbers can be exactly represented) gives 999999995904, i.e., the second operand is absorbed. The maximum error committed by representing a real number with a floatingpoint number under some rounding mode can be expressed in terms of the function ulp : R → R [START_REF] Muller | On the definition of ulp(x)[END_REF]. Its value on 1.0 is about 10 -7 for the single precision format.

The chop and near rounding modes are symmetric, i.e., the value after rounding does not depend on the sign: for each x ∈ R, •(x) = -•(-x).

B. Notation

R denotes the set of real numbers while F p,emax denotes an idealized set of binary floating-point numbers, defined from a given IEEE 754 format: this excludes subnormals and NaNs, but includes -∞, +∞ and zeroes. This restriction allows to considerably simplify the presentation (e.g., avoiding all technical details concerning subnormals); yet, everything can be generalized to any IEEE 754 binary floating-point format [START_REF] Bagnara | Filtering floating-point constraints by maximum ULP[END_REF]. The exposition is also much simplified by allowing e max to be ∞, i.e., by considering an idealized set of floats where the exponent is unbounded. F + p,emax denotes the "positive" subset of F p,emax , i.e., with s = 0. When the format is clear from the context, a real decimal constant (such as 10 12 ) denotes the corresponding round-tonearest, tail-to-even floating-point value (i.e., 999999995904 for 10 12 ). 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. Of course, we have

f max + = +∞ and (-f max ) -= -∞.
Binary arithmetic operations over the floats will be denoted by ⊕, ⊖, ⊗ and ⊘, corresponding to +, -, • and / over the reals, respectively. According to IEEE 754, they are defined with the rounding operator

• by x ⊕ y = •(x + y), x ⊖ y = •(x -y), x ⊗ y = •(x • y) and x ⊘ y = •(x/y).
As IEEE 754 floating-point numbers are closed by negation, we will denote the negation of x ∈ F p,emax simply by -x. The symbol ⊙ denotes any of ⊕, ⊖, ⊗ or ⊘. A floating-point variable x is associated to an interval of possible floatingpoint values; we will write x ∈ [x, x], where x and x denote the smallest and greatest value of the the interval, x ≤ x and either x = +0 or x = -0.

III. BACKGROUND ON CONSTRAINT SOLVING OVER FLOATING-POINT VARIABLES A. Interval-based Consistency on Arithmetic Constraints

Program analysis usually starts with the generation of an intermediate code representation in a form called three-address code (TAC). In this form, complex arithmetic expressions and assignments are decomposed into sequences of assignment instructions of the form result := operand 1 operator operand 2 . A further refinement consists in the computation of the static single assignment form (SSA) whereby, labeling each assigned variable with a fresh name, assignments can be considered as if they were equality constraints. For example, the TAC form of the floating-point assignment z

:= z * z + z is t := z * z; z := t + z, which in SSA form becomes t 1 := z 1 * z 1 ; z 2 := t 1 + z 1 ,
which, in turn, can be regarded as the conjunction of the constraints

t 1 = z 1 ⊗ z 1 and z 2 = t 1 ⊕ z 1 .
In an interval-based consistency approach to constraint solving over the floats, constraints are used to iteratively refine the intervals associated to each variable. A projection is a function that, given a constraint and the intervals associated to two of the variables occurring in it, computes a possibly refined interval for the third variable (the projection is said to be over the third variable). Taking z 2 = t 1 ⊕ z 1 as an example, the projection over z 2 is called direct projection (it goes in the same sense of the TAC assignment it comes from), while the projections over t 1 and z 1 are called indirect projections. Non-optimal projections for the four arithmetic operations can be found in [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF], [START_REF] Botella | Symbolic execution of floating-point computations[END_REF]. 4

B. The Marre-Michel Property

In [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF], Marre and Michel published an idea to improve the filtering of the addition/subtraction projectors. This is based on a property of the distribution 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 e x , if ∆ = x + -x, then for y of exponent e x + 1 we have This property, which can be generalized to subnormals, can intuitively be explained on Figure 1 as follows. Let z ∈ F p,∞ be a strictly positive constant such that z = x ⊖ y, where x, y ∈ F p,∞ are unknown. The Marre-Michel property says that y cannot be greater than α. In fact, α is carefully positioned so that α ++ -α + = 2(α + -α), e α + 1 = e β and z = β -α; if we take y = α + we need x > β if we want z = x -y; however, the smallest element of F p,∞ that is greater than β, β + , is 2∆ away from β, i.e., too much. Going further with y does not help: if we take y ≥ α + , then y -α is an odd multiple of ∆ (one ∆ step from α to α + , all the subsequent steps being even multiples of ∆), whereas for each x ≥ β, x -β is an even multiple of

y + -y = 2∆. Proposition 3.1: [20, Proposition 1] Let z ∈ F p,∞ be such that 0 < z < +∞; let also z = 1.b 2 • • • b i k 0 • • • 0 × 2 ez , with b i = 1; α = 1.1 • • • 1 × 2 ez+k , with k = p -i; β = α ⊕ z. Then, for each x, y ∈ F p,∞ , z = x ⊖ y implies x ≤ β and y ≤ α. Moreover, β ⊖ α = β -α = z. z α α + α ++ β β + z ∆ 2∆ 2∆ 0
∆. Hence, if y > α, z -(x-y) ≥ ∆ = 2 ez+1-i . However, since k = p -1, z + -z = z -z -= 2 ez+1-p ≤ ∆. The last inequality, which holds because p ≥ i, implies z = x ⊖ y.
A similar reasoning allows to see that x cannot be greater than β independently from the value of y.

In order to improve the filtering of the addition/subtraction projectors, in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF], Marre and Michel presented an algorithm to maximize the values of α and β over an interval. That algorithm and the main ideas behind the work presented in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] will be revisited and discussed in detail in Section IV-C.

IV. FILTERING BY MAXIMUM ULP

This section reformulates the Marre-Michel property so as to generalize it to multiplication and division. The filtering algorithms that result from this generalization are collectively called filtering by maximum ULP.

A. Upper Bound

For each floating-point operation ⊙ ∈ {⊕, ⊖, ⊗, ⊘}, we will define the sets F ⊙ ⊆ F p,emax and F ⊙ ⊆ F + p,∞ . Then we will define a function δ ⊙ :

F ⊙ → F ⊙ satisfying the following properties, for each z ∈ F ⊙ \ {-0, +0, -∞}: ∃y ∈ F ⊙ . δ ⊙ (z) ⊙ y = z; (1) ∀z ′ ∈ F ⊙ : z ′ > δ ⊙ (z) =⇒ ∄y ∈ F ⊙ . z ′ ⊙ y = z. (2)
In words, δ ⊙ (z) is the greatest float in F ⊙ that can be the left operand of ⊙ to obtain z. Remark that we may have F ⊙ F p,emax : properties (1) and (2) refer to an idealized set of floating-point numbers with unbounded exponents.

Since we are interested in finding the upper bound of δ ⊙ (z) for z ∈ [z, z], we need the following Proposition 4.1: Let w, v 1 , . . . , v n ∈ F ⊙ \{-0, +0, -∞} be such that, for each i = 1, . . . , n, δ ⊙ (w) ≥ δ ⊙ (v i ). Then, for each i = 1, . . . , n and w ′ ∈ F ⊙ such that w ′ > δ ⊙ (w), there does not exist a float y ∈ F ⊙ such that w ′ ⊙ y = v i .

Let z = x⊙y be a constraint where

-0, +0, -∞ / ∈ [z, z] and let w ∈ [z, z] be such that δ ⊙ (w) ≥ δ ⊙ (v) for each v ∈ [z, z]: then no element of x that is greater than δ ⊙ (w)
can participate to a solution of the constraint.

Dually, in order to refine the upper bound of y subject to the constraint z = x ⊙ y, it is possible to define a function δ ′ ⊙ that satisfies properties similar to (1) and [START_REF] Weyuker | On testing non-testable programs[END_REF]. In this paper we will focus on bounds for x only. Note, though, that when ⊙ is commutative (i.e., ⊕ or ⊗), δ ⊙ = δ ′ ⊙ .

B. Lower bound

For computing the lower bound, we will introduce functions δ ⊙ :

F ⊙ → F ⊙ for each z ∈ F ⊙ \ {-0, +0, +∞}: ∃y ∈ F ⊙ . δ ⊙ (z) ⊙ y = z; (3) ∀z ′ ∈ F ⊙ : z ′ < δ ⊙ (z) =⇒ ∄y ∈ F ⊙ . z ′ ⊙ y = z. (4)
These properties entail a result similar to Proposition 4.1: given the constraint z = x ⊙ y where -0, +0, +∞ / ∈ [z, z] and w ∈ [z, z] such that δ ⊙ (w) ≤ δ ⊙ (v) for each v ∈ [z, z], the float δ ⊙ (w) is a possibly refined lower bound for x.

C. Filtering by Maximum ULP on Addition/Subtraction

In this section we introduce the functions δ ⊕ and δ ⊕ . The functions δ ⊖ and δ ⊖ can be deduced by symmetry, as explained in Section III-B and [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF]. Using the Michel and Marre property (Proposition 3.1) we formally define the function δ ⊕ as follows.

Definition 4.1: Let us define

F ⊕ = F p,emax , F ⊕ = F + p,∞ , and let z ∈ F ⊕ be such that |z| = 1.b 2 • • • b i 0 • • • 0 × 2 ez , with b i = 1. Similarly to Proposition 3.1, let k = p -i, α = 1.1 • • • 1 × 2 ez+k and β = α ⊕ z. Then δ ⊕ : F ⊕ → F ⊕ is defined, for each z ∈ F ⊕ , as follows: δ ⊕ (z) =          +∞, if z = -∞ or z = +∞; α, if -∞ < z < 0; +0, if z = -0 or z = +0; β, if 0 < z < +∞.
Theorem 4.1: δ ⊕ is well-defined and satisfies (1) and ( 2). The function δ ⊕ : F ⊕ → F ⊕ is defined dually: for each z ∈ F ⊕ \ {-0, +0, +∞}, δ ⊕ (z) = -δ ⊕ (-z). It is easy to see that properties (1) and (2) of δ ⊕ entail properties (3) and (4) of δ ⊕ .

We now need algorithms to maximize δ ⊕ and minimize δ ⊕ over an interval of floating-point values. Since the two problems are dual to each other, we will focus on the maximization of δ ⊕ . As δ ⊕ is not monotonic, a nontrivial analysis of its range over an interval is required. When the interval contains only finite, nonzero and positive (resp., negative) values, the range of δ ⊕ has a simple shape. We are thus brought to consider an interval [z, z] such that z / ∈ {-∞, -0, +0} and z / ∈ {-0, +0, +∞} have the same sign. We will now revisit and correct the algorithm proposed by Michel and Marre in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] to maximize δ ⊕ over [z, z].

The idea presented in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] is the following. When dealing with an interval [z, z] with z > 0, α (and thus β and our δ ⊕ ) grows with the exponent and the number of successive 0 bits to the right of the mantissa, i.e., k in Proposition 3.1 and in Definition 4.1. Thus, maximizing these two criteria allows to maximize α over the interval.

Definition 4.2: Let z be a variable over F p,emax such that z / ∈ {-∞, -0, +0} and z / ∈ {-0, +0, +∞} have the same sign and z < z. Then µ ⊕ (z) ∈ F p,emax is given by As we have already pointed out, the algorithm of Definition 4.2 is very similar to the algorithm presented in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF].

1.0 • • • 0 × 2 ez , if e z = e z ; otherwise we define µ ⊕ (z) = 1.b 2 • • • b i a0 • • • × 2 ez , where, for some b i+1 = b ′ i+1 , |z| = 1.b 2 • • • b i b i+1 • • • × 2 ez ; |z| = 1.b 2 • • • b i b ′ i+1 • • • × 2 ez ; a = 0, if 1.b 2 • • • b i 0 • • • 0 × 2 ez = |z|; 1, otherwise.
There is an importance difference, though: in the case when

z = 1.b 2 • • • b i b i+1 0 • • • 0 × 2 ez , z = 1.b 2 • • • b i b ′ i+1 • • • × 2 ez and b j = 1,
for some j ≤ i. In this case, the algorithm of [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF]

returns 1.b 2 • • • b i 10 • • • 0 × 2 ez . Note, however, that the value that maximizes α is z, which is different from 1.b 2 • • • b i 10 • • • 0 × 2 ez .
Definition 4.2 can be extended to intervals that include subnormals [START_REF] Bagnara | Filtering floating-point constraints by maximum ULP[END_REF], but not to intervals containing zeroes. So, when z's interval contains zeroes, only the classical filtering is applied. For efficiency reasons, filtering by maximum ULP is only applied when δ ⊕ µ ⊕ (z) ≤ f max so as to avoid the use of wider floating-point formats. This example shows that filtering by maximum ULP can be stronger than classical interval-consistency based filtering. However, there are trivial examples that show the opposite phenomenon so that classical and maximum ULP are orthogonal: both should be applied, through interval intersection, in order to obtain optimal results.

D. Filtering by Maximum ULP on Multiplication

Let z ∈ F p,emax be a strictly positive constant such that z = x ⊗ y, where x, y ∈ F p,emax are unknown. As for Property 3.1, there exists a greatest float x m ∈ F p,emax such that there exists y ∈ F p,emax satisfying z = x m ⊗ y. More precisely, x m is the float such that z = x m ⊗ f nor min . Such a float always exists because multiplication by f nor min = 2 emin is equivalent to an exponent shifting. Now let us consider x ′ ∈ F p,emax such that x ′ > x m . By monotonicity of ⊗, z < x ′ ⊗ f nor min and there is no other float y = f nor min such that z = x ′ ⊗ y. In fact, by monotonicity, such float y should be smaller than f nor min . On the other hand, y must be greater than +0 for otherwise x ′ ⊗ y would not be strictly positive. However, for no y ∈ F p,emax we have +0 < y < f nor min . Therefore, the value x m such that z = x m ⊗ f nor min is the greatest value for x that can satisfy z = x m ⊗ y for some y.

Definition 4.3:

F ⊗ = z ∈ F p,emax f nor min • |z| ≤ f max and F ⊗ = F +
p,emax are the domain and codomain of δ ⊗ : F ⊗ → F ⊗ , defined for each z ∈ F ⊗ as follows:

δ ⊗ (z) = |z| • 2 -emin
Theorem 4.3: Function δ ⊗ is well-defined and satisfies (1) and ( 2).

The function δ ⊗ is simply defined as δ ⊗ = -δ ⊗ (z). Moreover, Definition 4.3 and Theorem 4.3 can be extended to intervals that include subnormals replacing all occurrences of f nor min by f min [START_REF] Bagnara | Filtering floating-point constraints by maximum ULP[END_REF]. The value M ∈ [z, z] that maximizes δ ⊗ is the one with the greatest absolute value, i.e., M = max |z|, |z|}. Since δ ⊗ is defined asδ ⊗ (z), the value that minimizes δ ⊗ is again M . Hence, if [z, z] does not contain zeroes, δ ⊗ (M ) (resp., δ ⊗ (M )) is an upper bound (resp., a lower bound) of x w.r.t. the constraint z = x ⊗ y. The restriction to intervals not containing zeroes is justified by the fact that, e.g., if z = 0 then z = x ⊗ y holds with x = f max and y = 0, hence no useful filtering can be applied to x.

As the product is commutative, the function of Definition 4.3 can be used for filtering y as well. Note that this filtering can only be applied when max |z|, |z| ∈ F ⊗ .

Example 4.2: Consider the IEEE 754 single-precision

constraint z = x ⊗ y with z ∈ [1.0 × 2 -50 , 1.0 × 2 -30 ] and x, y ∈ [-∞, +∞]. We have δ ⊗ (1.0 × 2 -30 ) = 1.0 × 2 -30 • 2 -(-126) = 1.0 × 2 96 , δ ⊗ (1.0 × 2 -30 ) = -1.0 × 2 -30 • 2 -(-126) = -1.0 × 2 96
, so, while classical filtering does not prune the intervals for x and y, filtering by maximum ULP yields the refined intervals

x, y ∈ [-1.0 • • • 0 × 2 96 , 1.0 • • • 0 × 2 96 ].

E. Filtering by Maximum ULP on Division

On the F p,emax domain, a role similar to the one of f nor min in the definition of filtering by ULP max on multiplication is played by f max in the definition of filtering by ULP max on division.

Definition 4.4: Let

F ⊘ = z ∈ F p,emax |z| ⊗ f max ≤ f max and F ⊘ = F + p,emax . Then δ ⊘ : F ⊘ → F ⊘ is defined, for each z ∈ F ⊘ , by δ ⊘ (z) = |z| ⊗ f max .
Theorem 4.4: Function δ ⊘ is well-defined and satisfies ( 1) and ( 2).

A similar result can be obtained for intervals that include subnormals by a suitable modification of Definition 4.4 [START_REF] Bagnara | Filtering floating-point constraints by maximum ULP[END_REF].

The function δ ⊘ is simply defined, for each z ∈ F ⊘ , by

δ ⊘ = -δ ⊘ (z).
The value M ∈ [z, z] that maximizes δ ⊘ is the one that has the greatest absolute value, i.e., M = max |z|, |z| . Since δ ⊘ is defined asδ ⊘ (z), M is also the value that minimizes δ ⊘ . Hence, if [z, z] does not contain zeroes, δ ⊘ (M ) (resp., δ ⊘ (M )) is an upper bound (resp. a lower bound) of x w.r.t. the constraint z = x ⊘ y. Once again, the restriction to intervals not containing zeroes is justified by the fact that, e.g., if z = 0 then z = x ⊘ y holds with x = f max and y = ∞, hence, also in this case, no useful filtering can be applied to x. Note that this filtering can only be applied when max |z|, |z| ∈ F ⊘ .

δ ⊘ (1.0 × 2 -110 ) = 1.0 × 2 -110 • 1.1 • • • 1 × 2 127 = 1.1 • • • 1 × 2 17 , δ ⊘ (1.0 × 2 -110 ) = -1.0 × 2 -110 • 1.1 • • • 1 × 2 127 = -1.1 • • • 1 × 2 17 .

Again, filtering by maximum ULP improves upon classical filtering with

x ∈ [-1.1 . . . 1 × 2 17 , 1.1 . . . 1 × 2 17 ].

F. Synthesis

Table I provides a compact presentation of filtering by maximum ULP under the assumption e min = 1 -e max , where the required functions can be summarized as follows:

δ ⊕ (z) = β, if 0 < z < +∞, α, if -∞ < z < 0; δ ⊕ (z) = -δ ⊕ (-z); δ ⊗ (z) = |z| • 2 -emin ; δ ⊗ (z) = -δ ⊗ (z); δ ⊘ (z) = |z| ⊗ f max ; δ ⊘ (z) = -δ ⊘ (z).

V. IMPLEMENTATION IN FPSE

A. FPSE FPSE [START_REF] Botella | Symbolic execution of floating-point computations[END_REF] is a constraint solver based on interval consistency filtering dedicated to the analysis of IEEE 754 floatingpoint computations coming from C programs. The tool takes a path condition as input, which is a quantifier-free conjunction of constraints extracted from a path of a C function. Constraints hold over the input variables of the program, including global variables, as well as temporary variables introduced by classical compiler code transformations. For a given path condition, FPSE can either return the first solution found or show that there is no solution. In the former case, the result can be interpreted as a test data that activates the selected path; in the latter case, infeasibility of the path is proved. Of course, solving the constraints in reasonable time is not always possible since the search space can be huge.

The constraints in FPSE are based on expressions built over ⊕, ⊖, ⊗, ⊘ and the relations =, =, <, ≤. Interval constraints (e.g., x ∈ [a, b]) are allowed as well as IEEE 754 type casting constraints, namely float-to-double, double-tofloat, long-to-double and double-to-long. FPSE works under some hypotheses that are now summarized. The tool deals only with the near tail-to-even rounding mode, which is used by default in almost all C implementations and is the most difficult to handle in constraint solving over floatingpoint variables [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF]. To model floating-point computations, special attention is paid to conform to the actual execution of programs. In order to capture the semantics of the program, it is of course necessary to respect the precedence of expression operators as specified by the C language as well as the evaluation order realized by the language implementation at hand. FPSE respects the shape of expressions as represented in the abstract syntax tree of the program without any rearrangement or simplification. The order in which operands are evaluated by a C implementation can be matched by using a preprocessor like CIL [START_REF] Necula | CIL: Intermediate language and tools for analysis and transformation of C programs[END_REF].

Any symbolic expression is decomposed into a sequence of TAC assignments where fresh temporary variables are introduced bearing in mind that the order of evaluation must be preserved. 5 This decomposition requires that intermediate results of an operation conform to the type of storage of its operands. 6 Constraint solving is implemented by using interval consistency combined with search heuristics. Several heuristics with static and dynamic choice of variable have been considered. FPSE is implemented with about 10 KLOC of SICStus Prolog (for the high-level constraint-solving machinery) and C (for the projection functions).

Table I FILTERING BY MAXIMUM ULP

SYNOPSIS Constraint x ⊆ • y ⊆ • Condition(s) z = x ⊕ y, 0 < z ≤ fmax [ δ ⊕ (ζ ), δ ⊕ (ζ )] [ δ ⊕ (ζ ), δ ⊕ (ζ )] ζ = µ ⊕ (z), -fmax ≤ δ ⊕ (ζ), δ ⊕ (ζ) ≤ fmax z = x ⊕ y, -fmax ≤ z < 0 [-δ ⊕ (ζ ′ ), -δ ⊕ (ζ ′ )] [-δ ⊕ (ζ ′ ), -δ ⊕ (ζ ′ )] ζ ′ = µ ⊕ (-z), -fmax ≤ δ ⊕ (ζ ′ ), δ ⊕ (ζ ′ ) ≤ fmax z = x ⊖ y, 0 < z ≤ fmax [ δ ⊕ (ζ ), δ ⊕ (ζ )] [-δ ⊕ (ζ ), -δ ⊕ (ζ )] ζ = µ ⊕ (z), -fmax ≤ δ ⊕ (ζ), δ ⊕ (ζ) ≤ fmax z = x ⊖ y, -fmax ≤ z < 0 [-δ ⊕ (ζ ′ ), -δ ⊕ (ζ ′ )] [ δ ⊕ (ζ ′ ), δ ⊕ (ζ ′ )] ζ ′ = µ ⊕ (-z), -fmax ≤ δ ⊕ (ζ ′ ), δ ⊕ (ζ ′ ) ≤ fmax z = x ⊗ y, 0 < |z| ≤ 2(2 -2 1-p ) [ δ ⊗ (m), δ ⊗ (m)] [ δ ⊗ (m), δ ⊗ (m)] m = max |z|, |z| , z = x ⊘ y, 0 < |z| ≤ 1 [ δ ⊘ (m), δ ⊘ (m)] m = max |z|, |z|

B. Relative ǫ

Slow convergence phenomena typically arise when few values are continuously pruned in a constraint propagation cycle. For example, the constraint x < y ∧ y ≤ x causes a slow convergence phenomenon over floating-point (and, for that matter, also integer) variables. Each time a projection function is woken, a single float is pruned from the domain of x and y. Unsatisfiability is ultimately proved, but not in reasonable time. To avoid slow convergence phenomena, we implemented a procedure that stops the filtering under a given threshold called relative ǫ. For a given floating-point variable, if a filtering function does not reduce its domain of more than ǫ%, then we withdraw that filtering function and do not prune the domain. An important difference with respect to interval-propagation-based constraint solvers over continuous domains is that we differentiate the treatment of direct and indirect projection functions. For direct projection functions this threshold is positioned at 0%, while for indirect projection functions it is positioned at 10% (a value that was experimentally determined to provide a good compromise). The idea is to benefit from the structure of the problem. When all the input variables of the program under test are instantiated, direct projection functions are sufficient to get a solution: thus applying them unconditionally is advantageous. In contrast, when a local or an output variable is instantiated first, this may lead to a propagation cycle that has to be cut very early: we discovered that cutting it using the relative ǫ on indirect projection functions is very effective as every potential propagation cycle involves at least one indirect projection function. For example, in x < y∧y ≤ x only indirect projection functions are involved on both variables, and propagation cycles do not occur. Of course, in this pathological case there is no propagation cycle but the system is still partially consistent at the end of the initial filtering: hopefully other constraints will allow to prune many values from the domains of x and y so that enumeration will not have to try all the values to prove unsatisfiability.

VI. EXPERIMENTAL EVALUATION

The aim of our experimental study was to evaluate filtering by maximum ULP (in brief ULP Max) and to determine whether it is an effective, practical property for solving constraints over the floats with an acceptable overhead. For presentation of the results in this paper we selected two C functions performing intensive floating-point computations. The first one is a small C function that computes a root of a polynomial equation within a given range: dichotomic() in Figure 2. Its computations are dominated by singleprecision floating-point computations. The second program is a real-world program embedded on unmanned airplanes to avoid fly-to-fly collision.

We implemented several search heuristics with static and dynamic variable orderings. For the choice of values, we implemented a domain-splitting strategy adapted to floatingpoint variables that proved to be very effective. For a given variable, our strategy selects first the floating-point value v that separates the domain of x in two equivalently-sized sub-domains, i.e., two domains containing the same number of floats; then it considers x = v, x < v and x > v by successively backtracking on these choices.

For the dichotomic() function we selected at random a path that iterates 12 times in the loop and considered all its path prefixes from iteration 1 to 12. For each path, we used FPSE to automatically generate a test input (an instantiation of all the input variables) that covers the path. We considered two versions of FPSE: a version that implements ULP Max filtering as defined in this paper, and a version without that. We measured: the number of elementary (ternary) constraints on the path (NbC); the number of uninstantiated floating-point variables on the path and the number of variables involved in the solution path of the search tree (NbV); the number of times ULP Max filtering takes place, globally and on the solution path (NbE); the number of floats pruned by ULP Max, in millions, globally and on the solution path (NbD); the percentage of domains pruned by ULP Max over all the variables involved in the solution path (%); the CPU time for generating test inputs with the standard version of FPSE (w/o) and with FPSE augmented For CPU time, we took the average of 10 runs of the same test input generation process. All results were computed on a system equipped with an Intel Core 2 Duo 3.00 GHz and running Linux 2.6 with 4 GB of RAM.

The results for dichotomic(), are presented in Table II. These show that FPSE with ULP Max is effective enough to solve at least 3 constraint systems that the standard version cannot solve. Moreover, ULP Max is able to prune the domains of floating-variables in all the cases and CPU time gains are due to this extra pruning power. On 6 paths, FPSE with ULP Max was able to prune between 8% and 20% of the variable domains.

The second example is a real-world example extracted from a critical embedded system. The function tcas_periodic_task_1Hz(), an excerpt of which is presented in Figure 3, is the core of a TCAS system (Traffic Collision Avoidance System) embedded into unmanned aircrafts. 7 The system receives the speed and direction of other aircrafts and, based on floating-point computations, it modifies the speed and direction of the host aircraft in order to avoid collisions. This program is interesting because determining the feasibility of its paths is hard and requires precise reasoning on non-linear floating-point computations.

For tcas_periodic_task_1Hz(), we selected all the possible paths up to 5 iterations of the main loop. This corresponds to about 130 paths among which 51% were shown to be infeasible by FPSE, regardless of whether ULP Max was used or not. For the remaining 49% of feasible paths, ULP Max has effects on 27 paths. We generated test inputs for all these paths with both versions of FPSE (i.e., w and w/o ULP Max). The results, which are reported in Table III, show that ULP Max always prunes the search space by more than 40 millions of single-precision floatingpoint values. In other words, ULP Max effectively prunes the search space in most cases. However, this extra-pruning does not always result in the overall speedup of the test input generation process: with reference to Table III, when the speedup is below 1 the overhead of computing ULP Max is not compensated by the gains it offers. On the other hand, it must be observed that those cases are not too frequent (9 cases out of 27 in this example) and the implementation of ULP Max in FPSE has much room for improvement.

VII. DISCUSSION

In this paper we brought to completion the plan anticipated in [START_REF] Carlier | Filtering by ULP maximum[END_REF]. This is part of a long-term research effort concerning the correct, precise and efficient handling of floating-point constraints [START_REF] Michel | Solving constraints over floating-point numbers[END_REF]- [START_REF] Botella | Symbolic execution of floating-point computations[END_REF], [START_REF] Belaid | Boosting local consistency algorithms over floating-point numbers[END_REF]- [START_REF] Carlier | Filtering by ULP maximum[END_REF].

Other authors have considered using search-based test data generation with a specific notion of distance in their fitness function [START_REF] Lakhotia | AUSTIN: A tool for search based software testing for the C language and its evaluation on deployed automotive systems[END_REF], [START_REF] Lakhotia | FloPSy: Search-based floating point constraint solving for symbolic execution[END_REF]. For instance, search-based tools like AUSTIN and FloPSy can generate a test input for a specific path by evaluating the path covered by some current input with respect to a targeted path in the program. However, they cannot solve the constraints of path conditions, since: 1) they cannot determine unsatisfiability when the path is infeasible, and 2) they can fail to find a test input while the set of constraints is satisfiable. 8 Recently, Borges et al. [START_REF] Borges | Symbolic execution with interval solving and meta-heuristic search[END_REF] combined a search-based test data generation engine with the RealPaver interval constraint solver, which is well-known in the Constraint Programming community. Even though FPSE and RealPaver are based on similar principles, their treatment of intervals is completely different. While FPSE preserves the solutions over the floats, RealPaver preserves the solutions over the reals 8 The floating-point intensive programs shown in the previous section seem to be outside the reach of the search-based tool AUSTIN [START_REF] Lakhotia | AUSTIN: A tool for search based software testing for the C language and its evaluation on deployed automotive systems[END_REF]: for example, AUSTIN seems to die pseudo-randomly on dichotomic(); both for dichotomic()and tcas_periodic_task_1Hz() we were not able to produce any test input data within 2 days of CPU time. by making the appropriate choices in the rounding modes used for computing the interval bounds. [START_REF] Botella | Symbolic execution of floating-point computations[END_REF] contains small examples showing that an interval constraint solver over the reals can miss floating-point solutions to constraints over floating-point variables. However, as RealPaver can treat transcendental functions with high precision, the approach followed in [START_REF] Borges | Symbolic execution with interval solving and meta-heuristic search[END_REF] allows the generation of floating-point inputs for programs that use such functions in a nontrivial way, something that is outside the scope of this paper.

VIII. CONCLUSION

This paper concerns constraint solving over floating-point numbers and its application to automatic test data generation. Interval-based consistency techniques are very effective for the solution of such numerical constraints, provided precise and efficient filtering algorithms are available. We reformulated and corrected the filtering algorithm proposed by Marre and Michel in [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] for addition and subtraction. We proposed a uniform framework that generalizes the property identified by Marre and Michel to the case of multiplication and division. The main ideas of this article were roughly sketched in [START_REF] Carlier | Filtering by ULP maximum[END_REF]: in this paper they have been revised, corrected and extended. The new filtering algorithms have been implemented in the FPSE system and the experimental evaluation show that they definitely improve the state-ofthe-art of automatic test data generation for floating-point programs. Future work includes the exploration of other properties based on linearization of floating-point computations, such as those proposed in [START_REF] Belaid | Boosting local consistency algorithms over floating-point numbers[END_REF].

Figure 1 .

 1 Figure 1. An illustration of the Marre-Michel property

Theorem 4 . 2 :

 42 Let z be as in Definition 4.2. Then, for each z ∈ [z, z], δ ⊕ (z) ≤ δ ⊕ µ ⊕ (z) .

Example 4 . 1 :

 41 Consider z = x⊕y with z ∈ [1.0, 2.0], x ∈ [-1.0 × 2 50 , 1.0 × 2 50 ] and y ∈ [-1.0 × 2 30 , 1.0 × 2 30 ]. With classical filtering we obtain x, y ∈ [-1.0 × 2 30 , 1.0 × 2 30 ], whereas with filtering by maximum ULP we obtain the much tighter interval x, y ∈ [-1.1 • • • 1 × 2 24 , 1.0 × 2 25 ].

Example 4 . 3 :

 43 Consider the IEEE 754 single-precision constraint z = x ⊘ y with z ∈ [-1.0 × 2 -110 , -1.0 × 2 -121 ] and x, y ∈ [-∞, +∞]. We have

Figure 2 .

 2 Figure 2. The dichotomic() function

Figure 3 .

 3 Figure 3. An excerpt of tcas.c

Table II EXPERIMENTAL

 II RESULTS FOR dichotomic() (TIMEOUT = 30 MIN)

	#	NbC	NbV	Global results NbE NbD	NbV	On the solution path NbE NbD	%	ULP Max w/o w/	Speedup factor
	1	17	12	62	17,515	12	1	864	20.2	0.142	0.080	1.775
	2	31	22	3,948	484,128	22	0	0	0.00	12.326	3.536	3.486
	3	45	32	461	102,522	32	3	1,174	9.15	3.969	0.872	4.552
	4	59	42	544,377	9,208,097	42	0	0	0.00	timeout	847.778	∞
	5	73	52	510	158,716	52	5	1,895	8.86	2.370	1.506	1.574
	6	87	62	799	209,621	62	0	0	0.00	timeout	2.050	∞
	7	101	72	494	87,934	72	7	2,625	8.77	6.087	0.983	6.192
	8	115	82	timeout	timeout	timeout	timeout	timeout	0.00	timeout	timeout	∞
	9	129	92	258	83,166	92	9	3,338	8.67	2.352	0.978	2.405
	10	143	102	637	157,421	102	0	0	0.00	timeout	2.482	∞
	11	157	112	224	73,702	112	11	4,034	8.57	2.471	0.724	3.413
	12	171	122	635	153,318	122	0	0	0.00	4.924	2.642	1.864
	float f(float x) { return x * x -2.0F; }
	float dichotomic(float xL, float xR) {
	float xM = 1.0F;
	while ((xR -xL) > 0.0001F) {
	xM = (xR + xL) / 2.0F;
	if ((f(xL) * f(xM)) > -1.0F) { xL = xM; }
	else { xR = xM; } }
	return xM; }

This is the the well-known oracle problem[START_REF] Weyuker | On testing non-testable programs[END_REF].

E.g., for the constraint z = x ⊕ y we have z = x ⊕ y and z = x ⊕ y (direct), x = mid(z, z + ) ⊖ y and x = mid(z, z -) ⊖ y (1 st indirect), y = mid(z, z + ) ⊖ x and y = mid(z, z -) ⊖ x (2 nd indirect). Here, for finite x, y ∈ Fp,e max , we denote by mid(x, y) the number that is exactly halfway between x and y; note that either mid(x, y) ∈ Fp,e max or mid(x, y) ∈ F p+1,emax .

The introduction of temporary variables does not change the semantics of floating-point computations as long as it reflects the behavior of the compiler and of the floating-point unit.

This is not always true: e.g., on Intel's architectures based on the 387 floating-point coprocessor, registers have more precision than the IEEE

float and double types; this makes rounding unpredictable. Luckily, SSE instruction sets, which do not pose this problem, are superseding 387.

The complete source code is available at http://paparazzi.enac.fr

ACKNOWLEDGMENT

We are grateful to Abramo Bagnara (BUGSENG srl, Italy) for the many fruitful discussions we had on the subject of this paper, and to Paul Zimmermann (INRIA Lorraine, France) for the help he gave us proving a crucial result.


Sagar Sen

Arnaud Gotlieb
email: arnaud@simula.no

Testing a Data-intensive System with Generated Data Interactions The Norwegian Customs and Excise Case Study

Keywords: database systems, pairwise, T-wise, software product lines, database schema, entity-relationship diagram models, feature models, feature diagram, Alloy, testing

Testing data-intensive systems is paramount to increase our reliance on e-governance services. An incorrectly computed tax can have catastrophic consequences in terms of public image. Testers at Norwegian Customs and Excise reveal that faults occur from interactions between database features such as field values. Taxation rules, for example, are triggered due to an interaction between 10,000 items, 88 country groups, and 934 tax codes. There are about 12.9 trillion 3-wise interactions. Finding interactions to uncover specific faults is like finding a needle in a haystack. Can we surgically generate a test database for interactions that interest testers? We address this question with a methodology and tool Faktum to automatically populate a test database that covers all T-wise interactions for selected features. Faktum generates a constraint model of interactions in Alloy and solves it using a divide-and-combine strategy. Our experiments demonstrate scalability of our methodology and we project its industrial applications.

Introduction

Data-intensive software are increasingly prominent in driving global processes such as scientific/medical research, E-governance, and social networking. Large amounts of data is collected, processed, and stored by these systems in databases. For example, the Norwegian Customs and Excise department uses the TVINN system to processes about 30,000 declarations a day. TVINN stores validated transactional information such as declarations in a central database. It processes incoming declarations to verify their conformance to well-formedness rules, customs laws and regulations before accepting a declaration in the database. This scenario is prevalent in many data-intensive software systems dealing with transaction data which comprises of semi-structured/structured data in medium/high volume. Testing the data-intensive software systems in the industrial context is the subject of this paper. At Certus Centre, our objective is to identify and provide innovative solutions to thriving industrial problems. Certus is founded as a consortium with several industry and public administration partners and is primarily funded by the Research Council of Norway. We interact with industry partners in weekly meetings to understand and help solve their problems using state-of-the-art research findings and tools. We go by the motto-Industry is our lab!, owing to the fact that our scientific problems emerge from industry. This paper narrates the story of our collaboration with the Norwegian Customs and Excise department who bring to us the problem of testing the TVINN system.

The TVINN system at the Norwegian Customs and Excise department has been in operation since 1988. It is a massive database application that processes declarations sent as standard EDIFACT messages from customs officers and companies concerned with import/export in Norway. A sophisticated script called EMIL, written in the legacy language Sysdul, processes about 30,000 declarations/day to notify customs officers about the correctness of the declarations. The correctness is verified based on a large set of well-formedness rules and customs laws and regulations. For instance, if the net weight is less than half the gross weight the script identifies a problem and sends back a declaration for further evaluation by a company/customs officer. The customs laws change periodically with new governmental policies. For example, there has been a steep hike in taxes for import of cheese from Europe to protect Norwegian farmers. Customs laws also can be change for a short period of time by Customs officers who want to enforce a law on a series of suspicious trucks. Therefore, the dynamic change of rules in the TVINN system makes testing it constant effort. Moreover, a recent major development in the TVINN system is its migration from Sysdul to a system developed in Java. Will the new system function equivalent to the one currently under operation? Will the system regress? This is the fear that a public service such has to face. Current practice in testing involves domain experts manually creating test databases that aim to reveal faults and differences between different versions of TVINN. Talking to testing experts we understand that manual test creation is tedious and most importantly its hard to know where the tester's journey ends. There is a lack of the notion of test coverage. In this paper, we outline Certus center's efforts to address this particular challenge.

We present a methodology to automatically generate test data that cover all valid T-wise interactions between database features. These features correspond to selected field values in the specification of a database's input domain which is the database schema. We use the feature modelling formalism to help a test engineer select variation points in a database schema as a feature model. It also allows the tester to specify forbidden feature combinations and feature interdependencies. Using a feature model as input we apply our previously developed approach in [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] to generate a set of configurations that cover all valid T-wise interactions between database fields of choice. This configuration generation module is scalable to large feature models and satisfies all feature interdependencies in the configuration. The approach is based on transforming a feature model and T-wise combinations to a set of constraint systems in Alloy. These systems are solved concurrently and based on a divide-and-combine strategy. The solutions are concurrently generated as sets of configurations that are combined into a final set. We transform these configurations to insertion queries to populate a test database. However, certain fields are still placeholders for values in their appropriate domains. The fields that have little significance to a tester are updated using generated SQL update queries. The update queries either update values to maintain data integrity (domain and referential integrity) or random values for fields with less significance to the tester. The result is a complete set of queries to populate a test database covering all valid T-wise interactions between database features of interest. We implemented the methodology in a prototype tool called Faktum (fact in Norwegian).

We performed experiments to validate our methodology. We used a common database schema developed in collaboration with the Norwegian Customs and Excise as the input domain and specify a feature model to select data features concerning imports from Brazil, India, China, and the USA. We generated a test database using these inputs that covers all 2-wise/pairwise interactions between a total of 75 features such as currencies, country codes, country groups, tax fee codes, and declaration categories to name a few. We generated about 935 configurations covering all valid pairwise interactions (of 10804 all possible pairwise interaction) between the chosen features. These 935 configurations were transformed to a set of Structured Query Language (SQL) queries to populate a test database. The time complexity to detect the validity of a single pair (/tuple) is O(1). The average time to detect validity is about 12 ms. We generated 187 sets of configurations covering divided subsets of pairs. The generation of configurations is time and scope bounded. Every call to Alloy's SAT solver is bounded to an average of 400 ms and 6803 calls to the solver were necessary to obtain the final set. The large number of calls to Alloy's SAT solver is a trade-off to address scalability. We combined a total of 935 configurations from 187 sets that were solved concurrently in 48 minutes.

We summarize our contributions as follows:

Contribution 1: We propose a scalable methodology and a prototype tool Faktum to automatically generate databases that cover all T-wise interactions between data fields of testing interest. Contribution 2: We demonstrate through experiments on an industrial case study that our methodology is scalable and can be used surgically generate tests.

The paper is organized as follows. In Section 2, we present a detailed overview of the case study at the Norwegian Customs and Excise and the problems associated with the complexity of the testing task. In Section 3, we present a methodology based on automatic generation of test databases satisfying T-wise interactions to address the testing problem at the Norwegian Customs and Excise. We present results from a concrete experiment in Section 4. Section 5 discusses of related work, while we conclude with a summary of our experience in Section 6. We describe our industrial case study from the Norwegian Customs and Excise as illustrated in Figure 1(a). The system under study is the Tollvesenets startside for internettfortolling (TVINN). An overview of TVINN process flow is presented in Figure 1(a) and official description is available at https://fortolling.toll.no/Tvinn-Internett. Customs officers and industries associated with import/export create declarations at Norwegian ports of entry. These declarations are encapsulated in the EDIFACT standard for business communication. A declaration is encapsulated as an EDIFACT CUSDEC message. These messages are sent to TVINN's central server where they are processed by a sophisticated script called EMIL. EMIL parses EDIFACT messages and verifies them against well-formedness rules. It then verifies if the declared amount is accurately computed based on a statistical value for an item. These rules depend on numerous factors such as (a) 260 countries of origin divided in 88 country groups (b) over 160 currencies (c) taxes are computed based on about 900 tax code groups (d) an list of more than 10,000 items. A declaration can be categorized into 6 different categories based on EMIL's computation. The simplest categories being complete and reject. The response from TVINN is sent back as an EDIFACT CUSRES message to customs officer or industry. Rules in TVINN evolve on a regular basis (approx. every 6 months) depending on new governmental policies, sanctions, and change in political parties. TVINN is also affected with time-bounded rules created by customs officers. These time-bounded rules are ephemeral and exist for a short period of time. For instance, a customs officer decides to thoroughly check 20 trucks coming from a nation X in civil war. He/she will possibly create a rule to check all trucks from X for the next 3 hours. This rule is called mask control will disappear after 3 hours. These rules can change on an everyday basis without anticipation making TVINN a highly dynamic system.

Industrial Case Study: Norwegian Customs and Excise

TVINN is a complex and dynamic database application that processes about 30,000 declarations per day. Testing TVINN is a challenge since it evolves rapidly. Moreover, in 2012-2013 TVINN will undergo a migration from its native implementation Sysdul to Java. Will the new implementation in Java regress with respect to the old one? This is the question that haunts the Norwegian Customs and Excise at present and in the years to come.

Testing TVINN has been intuitively achieved by a small testing staff executing a subset of a large number of readily available records of real-world declarations. However, using these real-world declarations or records present four important problems: No Coverage Guarantee: Records obtained from real-world transactions such as customs declarations cover a realistic subset of the database's domain (set of all possible combination of values in fields and tables of a database). However, they often do not cover combinations of values that are very rare or exceptional. Very Large Set of Test Records: Accumulating information from real-world transactions can easily give rise to an ever-growing set of data records. Many of these records share similarities and hence are redundant for the purpose of testing. Cost-effective testing will require a selection of a minimal set of records that precisely captures testing intentions. A minimal set will also have modest time and space requirements for testing efforts such as nightly tests. Confidentiality: Governments/enterprises involving financial transactions or military data for instance have stringent confidentiality agreements with their clients. Therefore, its often impossible for them to outsource their testing efforts to external agencies who would use real-world records. Constantly Changing Rules: Records for testing often have a lifetime and need to be discarded. For instance, in the Norwegian Customs and Excise system changes when sanctions are imposed on countries or significant changes happen in currency exchange rates. Legacy transaction records may not be used anymore to test the evolved system.

Our objective at Certus Software V&V center, Simula Research Laboratory, is to present a solution to address these problems.

Methodology

In this section, we present a methodology to generate test databases that cover all T-wise interactions between data fields of interest. The overview of our methodology is shown in Figure 3 (a). In Section 3.1, we present foundational notions that will be used to describe our methodology. We present the different steps in our methodology in Section 3.2. In Section 3.3, we describe the implementation of the methodology in our tool Faktum.

Foundations

Database Schema The first input artifact in our methodology is a database schema. It specifies the input domain of a database. We briefly describe the wellknown concept of a database schema. More information on them can be found in a standard database textbook such as in [START_REF] Date | An Introduction to Database Systems[END_REF]. A database schema typically contains one or more tables. A table contains fields with a domain for each field. Typical examples for field types/domains are integer, float, double, string, and date. The value of each field must be in its domain hence maintaining domain integrity in a database. A table contains zero or more records which is a set of values for all its fields within their domain. A table may contain one or more fields that are referred to as primary keys. This means that each record is identified by its primary key. Table may refer to primary keys in other tables via foreign keys. The value of foreign keys must match the value of a primary key in another table. This is known as a referential integrity constraint. We refer to the combined concepts of referential integrity and domain integrity as data integrity. Records in a database must satisfy data integrity as specified by its a database schema. Databases can be queried using Structured Query Language (SQL) queries. These queries are both used to create and populate a database and query it for information presented as a table or a view. In this paper, we generate several hundred SQL INSERT and UPDATE queries to populate a test database (see Section 3.2 and Section 4 for more information on generated queries).

As a running example, we present a schema developed along with our industry partner, the Norwegian Customs and Excise, in Figure 1(b). The database schema consists of four tables and is created on a MySQL server. We describe the tables and some of the fields in them. The Customers table is used to store records of customers. A customer is identified by a CustomerID which is a primary key (indicated PK). A customer can make one or more declarations. These declarations are stored in the Declarations table that refer to a customer using a foreign key (indicated FK). A declaration can have one or more items that is stored in the Items table. Every item has an item code and a statistical value of its cost. There can be different types of taxes on an item which is stored in the Taxes table. The most common form of tax is the value added tax or VAT. Taxation rules are often expressed on the country group, tax fee code group (from the Taxes table) of import and the item code (from the Items table). The 10,000 items codes, 88 country groups, and 934 tax fee codes can potentially give rise 12.9 trillion 3-wise possible taxation laws. However, only 195,000 taxation laws are used in practice.

Feature Model of Database Variability Populating the database schema requires selection amongst a set of choices for its field values. The second input artifact to our methodology is a model of variation points or choices in the database schema. We use the feature modelling formalism to specify the variability in a database schema. The feature modelling formalism is described in detail in [START_REF] Batory | Feature models, grammars, and propositional formulas[END_REF]. Typically, a feature model is used to specify the different features in a software product line and their inter-dependencies. Inter-dependencies are constraints on the choice of features. Some features are mandatory, some are optional, some features require other features while some features are mutually exclusive (XOR) with respect to other features. Features can be abstract or concrete. Abstract feature help in classification and hierarchy while concrete features go into a final software configuration. A configuration of a feature model is a finite set of concrete features that satisfy feature model constraints. The notion of feature models is very popular and a general formalism to specify variability in software product lines and software artifacts in general. We use the feature modelling formalism in this paper to specify variability in a database schema.

In Figure 2, we present the feature model of database variability in the Norwegian Customs and Excise schema. The mandatory root feature specifies the database identifier which is TollCustoms. All child features of root in the second level specify identifiers for the different tables available in the database. The third level contains features for database fields in each of the tables. In the fourth level of model the features specify the different possible values for each field. The database, table, and field features in the first, second, and third level respectively are mandatory but abstract features. The values for fields are mutually exclusive which means only one value can be associated to a field. For instance, Declarations is table feature that has a Category feature. There are six different declaration categories. Only one of the values such as FU can be associated to the category field. We present the field values in a concise manner in Figure 2. We show the first possible variation in a field value while we only show the number of other possible possible values due to space limitations. A selection of field value features for all fields across all tables is what we call a record configuration. A record configuration specifies the exact values that will go into each field for a record in each and every table of the database. We use the term configuration interchangeably with record configuration. Multiple record configurations represents multiple records that go into the different tables of the database. For instance, [F U, U SD, BR, Import, ...] is part of a record configuration that shows values that will go in as a record in the Declarations table .   Combinatorial Interaction Testing Combinatorial interaction testing (CIT) is a recognized software testing technique introduced by Cohen [START_REF] Cohen | The aetg system: An approach to testing based on combinatorial design[END_REF], that tests all interactions between features, parameter values or in our case database field vales. An interaction can be seen as a tuple of software features. A widely cited NIST study of the fault databases of several real-life systems reports that all the faults in these databases are caused by no more than six factors [START_REF] Kuhn | An investigation of the applicability of design of experiments to software testing[END_REF]. If test features are modelled effectively, a T-wise testing can expose all the faults that involve no more than T features. However, pairwise or 2-way testing has been shown to be both time efficient and effective for most real case studies http://www.pairwise.org/results.asp. This motivated us to focus on CIT of all pairwise interactions of database field values of interest (see Section 4). T-wise testing requires that every T-wise interaction of between all features/database field values is present at least once in a set of record configurations. Generating record configurations that cover all T-wise interactions is a challenging task for very large feature models. There is a combinatorial explosion in the number of possible interactions with the increase in the number of features and the value of T. For instance, 10,000 items codes, 88 country groups, and 934 tax fee codes give rise to 192.5 million pairwise/2-wise and 12.9 trillion 3-wise interactions. In [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF], we present a scalable approach and tool to generate a set of configurations that cover all T-wise interactions between features. We use the approach in [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] to generate a set of record configurations that cover all T-wise interactions between database field values. Although, our approach makes the problem tractable it still depends on the computing resources (running in parallel) for a very large number of interactions. Therefore, in our experiments we go one step further and surgically select only a relevant subset of field values.

Methodology Description

We describe the methodology in four phases where each phase subsumes several sub-steps. An illustrative overview of the methodology is given in Figure 3 (a). Phase 1. Tester interaction: A tester provides three inputs (a) A database schema specification as described in Section 3.1. (b) Variability in a database as a feature model. The variations in database field values for TollCustoms is shown in Figure 2 (c) The value of T which represents the strength of interactions needed to be tested. For instance, when T = 2 we intend to test all 2-wise or pairwise interactions between database field values. If we consider all pairwise interactions between two features declaration category F U and currency code U SD we have 2 2 different interactions as shown in Figure 3(b). Similarly, we have 2 3 possible 3-wise interactions between a set of 3 features. In general, with n features we have 2 T × n T possible interactions minus those interactions that are forbidden by feature model inter-dependencies.

However, we notice that the INSERT query only contains values for a subset of all fields in record. The remaining fields have a NULL value. Therefore, in Step 6, we generate UPDATE SQL queries to fill in values for the remaining fields. The values of the remaining fields can be pseudo-randomly generated in their respective domains while satisfying data integrity constraints. For instance, a value for CustomerID in the declaration table must be an numeric string with at least 8 digits (domain integrity) and every customer id must be present in the Customers table (referential integrity). In our methodology, coverage of Twise interactions has priority. Therefore, we generate unique random values for all remaining fields. All foreign key field values are identical to their primary key values to ensure referential integrity. In Listing 1.2, we illustrate an UPDATE SQL query to complete the partial record created using INSERT query in Listing 1.1.

Implementation in Faktum

We implement the database synthesis tool Faktum for T-wise interaction testing in Java. The implementation is standalone and uses libraries for constraint solving such as Alloy. A prototypical implementation of the tool to generate configurations is available online https://sites.google.com/a/simula.no/dbtwise/.

Experiments

We perform an experiment to synthesize a test database for the Norwegian Customs and Excise case study and discuss the database synthesis results, scalability, and threats to validity of our approach.

Experimental Setup

We develop a test scenario to automatically synthesize a database that will test rules applying to Norwegian imports/exports from/to Brazil, India, China and USA. The feature model of database variability is specified in Figure 4. This is a subset of the full feature model shown earlier in Figure 2. This feature model specifies the variations of fields of interest in the database schema shown in Figure 1(b). This selection of features surgically pin-points a specific testing zone in the vast input domain. We use the feature model in Figure 4 as input to generate a set of configurations that cover all 2-wise/pairwise interactions between field values. The total number of interactions is the number of interactions between 37 concrete features minus the invalid interactions due to the XOR relation between field values in the same field. There are 4

× 37 2 -5 2 -4 2 -4 2 -2 × 2 2 -4 2 -5 2 -4 2 -9 2
= 2582 pairwise interactions that need to be covered by a set of record configurations.

Results of Discussion

We generate 935 record configurations covering all valid 2582 pairwise interactions of features in Figure 4. In the worst case, a naive non-optimal set of configurations will have the size 2582 with one interaction covered per configuration. These record configurations are transformed to SQL queries to populate a test database. The set of SQL queries are available in https://sites.google.com/a/simula.no/dbtwise/.

In this paper, we focus on understanding the scalability of our approach for the industrial case study. An important step in our methodology verifies if an interaction between database field values is indeed consistent with respect to its feature model. Every interaction is a tuple of field values that must be present in a configuration. However, not all tuples are valid with respect to the feature model in Figure 4. For instance, a tuple cannot contain two different categories for a declaration since all categories are mutually exclusive (XOR). Therefore, these forbidden 1) for a feature model of finite size. This is experimentally validated in Figure 5. We measure the time spent in the block of code for checking tuple validity. We observe that the average time in this block is about 12 ms when measured every 10 seconds. We use the non-intrusive perf4J library to perform the measurements.

Another important step is the generation of a set of configurations satisfying all valid tuples of feature interactions. We use the divide-and-combine approach to concurrently generate configurations satisfying all interaction tuples divided in exclusive subsets. Our methodology creates 187 concurrent sets of configurations satisfying the 2582 interaction tuples. The total number of configurations in the 187 sets is 935. We measure the number of calls to the constraint solver as shown in Figure 6. Creating the 935 configuration required about 6804 calls to Alloy's SAT solver. However, the average time spent in the solver was within about 450 ms due to the finite scope of the solver. The large number of calls to a constraint solver is a trade-off to achieve high scalability. Theoretically, we can imagine just one call to a constraint solver attempting to generate a set of configurations satisfying all 2582 interactions. However, in practice this is intractable and our divide-and-combine strategy address this exact issue.

Threats to Validity

Our experiments do not consider a certain number of factors that could affect our generalized outcomes on scalability. The scalability of our methodology may most likely be affected by the increase in the number and complexity of constraints in the feature model. However, we suggest surgically creating test models covering interactions between only relevant features for a tester. This will most often keep the satisfaction problem tractable. Another, factor we did not consider Fig. 6. Concurrent Configuration Generation is presence of numerical constraints between database fields that are generated randomly. For instance, the value of computed tax for an item is a function of the statistical value of an item. In this paper, these database fields were associated to random values without conforming to numerical constraints since in the general intention for testing was the interaction and not he numerical correctness of fields storing real numbered values. An extension to our tool Faktum will be to ensure that numerical constraints are satisfied between real-valued fields.

Related Work

The related work for this paper falls into two areas: (a) generation of test configurations covering T-wise interactions in a feature model (b) generation of synthetic data for databases. We will review the some relevant essential contributions in these areas.

Researchers have proposed several approaches for generating test configurations satisfying T-wise and notably the 2-wise/pairwise coverage criteria. The pioneering work in this area is the AETG approach [START_REF] Cohen | The aetg system: an approach to testing based on combinatorial design[END_REF] and its implementation to address highly-configurable software systems [START_REF] Cohen | Constructing interaction test suites for highlyconfigurable systems in the presence of constraints: A greedy approach[END_REF]. Based on a greedy algorithm, AETG generates N-wise covering arrays for a set of parameters and converts these arrays into a set of test cases. However, it cannot deal with constraints among the parameters, thus limiting its adoption to feature models with constraints. Oster et. al. [START_REF] Oster | Automated incremental pairwise testing of software product lines[END_REF] uses a greedy and ad-hoc algorithm based on the maximum number of valid pairs within each configuration while satisfying constraints. Recently, Johansen [START_REF] Johansen | Properties of realistic feature models make combinatorial testing of product lines feasible[END_REF] introduced the SPLCATool (Software Product Line Covering Array Tool) to generate test configurations from feature models. Similar to the above mentioned approaches, the tool uses greedy algorithm to enforce all pairs in a set of configurations. The tool is quite efficient and has been used on large feature models. The tool PACOGEN developed by Hervieu et. al. [START_REF] Hervieu | Pacogen: Automatic generation of pairwise test configurations from feature models[END_REF] goes a step further and generates a minimal set of configurations that covers all pairwise using a time-aware constraint solving procedure. However, for the problem of generating database records covering all T-wise interactions it was necessary to populate a test database with the generated configurations. Therefore, we used our previously developed approach described in Perrouin et.al. [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] to generate a set of configurations covering T-wise and satisfying feature model constraints. The approach in [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] is scalable as it based on a divide-and-combine strategy. According to our knowledge, this is the first time an interaction testing tool is used to populate a test database for data-intensive systems.

Important work to generate synthetic databases include [START_REF] Bruno | Flexible database generators[END_REF], [START_REF] Houkjaer | Simple and realistic data generation[END_REF], [START_REF] Lo | Generating databases for query workloads[END_REF], [START_REF] Mannila | Automatic generation of test data for relational queries[END_REF]. Database generation tools such as in [START_REF] Houkjaer | Simple and realistic data generation[END_REF], [START_REF]IBM: IBM DB2 test data generator[END_REF], [START_REF] Bruno | Flexible database generators[END_REF] allow users to specify the data distributions over attributes and intra-attribute correlations. Houkjaer et al. [START_REF] Houkjaer | Simple and realistic data generation[END_REF] use a graph model containing primary-foreign keys. This model is used to guide the data generation process. In all these works, the question of generating general T-wise interactions-covering data is not addressed and considering constraints among the fields is irrelevant. In that respect, the methodology proposed in this paper innovates also w.r.t. to database generation tools.

Conclusion

In this paper, we address the problem of testing TVINN, a data-intensive system in the Norwegian Customs and Excise department. TVINN processes all custom declaration coming in and out of Norway. We introduce a methodology implemented in a tool called Faktum, to populate a test database with configurations that cover all T-wise interactions between selected custom declaration database fields. In an experiment, we develop a variability model (feature model) and a database schema both developed in collaboration with the Norwegian Customs and Excise department. We use our tool to automatically generate 935 configurations covering all valid 2582 2-wise/pairwise interactions of features in less than 48 minutes. The scalability of our approach and the possibility to generate test data covering all T-wise interactions has given way to food for thought in our academia-industry partnership. For instance, we realize that random elements, in the test database, have a low degree of comprehensibility for testers at TVINN. Hence, we propose the use of Faktum for partial population of test databases. The remaining elements are manually completed by testers at Norwegian Customs and Excise. The partial test data gives testers the confidence of covering all valid T-wise interactions. Our future work with our partner will involve a large-scale empirical evaluation of test generation while carefully taking into account the clauses of non-disclosure agreements.

Fig. 1 .

 1 Fig. 1. (a) Norwegian Toll Customs Industrial Case Study (b) Common Database Schema at Norwegian Toll Customs in Crow-Foot Notation

Fig. 2 .

 2 Fig. 2. Feature Model of Database Variability

Listing 1 . 2 .

 12 UPDATE INTO D e c l a r a t i o n s ( CustomerID , Date , S e q u e n c e , V e r s i o n , Amount , FeeAmount , T r a n s p o r t C o s t , E xc h a n g e R a te ) VALUES ( ' 2 0 0 2 5 4 2 6 1 6 ' , ' 1965-3-29 ' , ' 1 ' , ' 1 ' , ' 2 9 8 2 . 4 9 0 2 4 5 ' , ' 1 UPDATE query to Complete a Declarations Record Phase 4. Population of Test Database: The last phase of our methodology is a straightforward population of a MySQL test database. The database is populated using INSERT queries followed by UPDATE queries (as shown in Step 8, Figure3(a)) generated in Phase 3 of our methodology. The test database generated in for the common database schema can be exported to more sophisticated industrial data-intensive systems in Step 9. This separation of database specifications allow us to work on a testing-specific subset of the entire database at Norwegian Customs and Excise. The full TVINN database contains several entities and fields that do not interest testers. The collaborative effort with the Norwegian Customs and Excise department helped us extract a testing specific subset with bi-directional exportability of records.

Fig. 4 .

 4 Fig. 4. Feature Model with Features Selected for Interaction Testing

Fig. 5 .

 5 Fig. 5. Checking Tuple Validity has O(1) Time Complexity

Acknowledgement

We thank the Norwegian Customs and Excise department for their trustful interactions with us. In particular, we would like to thank Atle Sander, Astrid Grime and Katrine Langset for their valuable inputs. We thank the Research Council of Norway for their generous support ; it would not have been impossible to setup such a close industry-academia collaboration for high impact software engineering research without it.

 2. This is achieved via steps enumerated 1-4 in Figure 3(a). These steps are derived from our previous work in [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF]. We briefly describe the approach from [START_REF] Perrouin | Automated and scalable t-wise test case generation strategies for software product lines[END_REF] here. In Step 1, we automatically transform a feature model to a constraint satisfaction problem A in the formal language Alloy [6] [START_REF] Jackson | Software Abstractions: Logic, Language, and Analysis[END_REF]. In Step 2, we generate Alloy predicates that encode tuples of features representing T-wise interactions between concrete field value features (shown in grey in Figure 2). We insert these predicates into the Alloy model A. In Step 3, we detect all T-wise Alloy predicates consistent with A and reject the others which are not accepted by the feature model specification.

In Step 4, we use a divide-and-combine strategy to generate sets of configurations that satisfy interaction tuples divided in subsets. We combine the sets to obtain a set of configurations that cover all T-wise interactions between features.