Claude-Pierre Jeannerod

Jean-Michel Muller

Paul Zimmermann

On various ways to split a floating-point number

We review several ways to split a floating-point number, that is, to decompose it into the exact sum of two floatingpoint numbers of smaller precision. All the methods considered here involve only a few IEEE floating-point operations, with rounding to nearest and including possibly the fused multiply-add (FMA). Applications range from the implementation of integer functions such as round and floor to the computation of suitable scaling factors aimed, for example, at avoiding spurious underflows and overflows when implementing functions such as the hypotenuse.

I. INTRODUCTION

Splitting a floating-point number (that is, separating a precision-p floating-point number x into two floating-point numbers x h and x of smaller significand size, and such that x = x h + x ) has several interesting applications in floatingpoint arithmetic. Among them (the examples below are in radix 2 arithmetic):

• it can be used to implement various integer functions on floating-point numbers: for example, functions such as round(x) (round x to a nearest integer) and x/2 k are special cases of splittings; • it can be used for accurately computing sums of several floating-point numbers. The main idea is to split these numbers so that the obtained values can be affected to containers (also called bins), such that all floating-point numbers in the same bin have significand of at most b bits and have their least significant bit of a fixed weight (that depends on the bin). This allows for error-free accumulation of all numbers of the bin provided that we do not add more than 2 p-b numbers to one bin. This method seems to appear first in a Matlab program in [START_REF] Zielke | Genaue Lösung linearer Gleichungssysteme[END_REF]. It is presented in detail and analyzed in [START_REF] Rump | Accurate floating-point summation, Part I: Faithful rounding[END_REF]; • it can be used for computing the error of a floatingpoint multiplication on architectures that do not offer an efficient FMA (fused multiply-add). The underlying idea is to split each of the operands into two p/2 -bit numbers, so the resulting numbers can be multiplied together without error. The first presentation of this method goes back to [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], and the splitting it uses is attributed to Veltkamp; • it can be used for computing a power of 2 close to |x|.

This is useful for scaling some calculations. Consider for instance the evaluation of √ a 2 + b 2 in floatingpoint arithmetic. The computation of a 2 or b 2 (or their sum) may underflow or overflow, resulting in a very inaccurate or infinite final result, even when the exact value of √ a 2 + b 2 is far from the underflow and overflow thresholds. A frequently suggested way of dealing with the problem is to divide both operands by the scaling factor max{|a|, |b|}, which gives new values a and b , to compute √ a 2 + b 2 , and to multiply the obtained value by the same scaling factor. Unfortunately, dividing and multiplying by the scaling factor is in general not exact. A way to avoid this is to choose as scaling factor an integer power of 2 close to max{|a|, |b|}.

In the first two examples, we wish to split numbers around a constant, that is, we want x h to be multiple of some constant (1 in the case of function round), and |x | to be less than this constant (or less than half this constant). We will call these splittings absolute splittings. In the last two examples, we wish to split numbers around a value that is somehow proportional to |x| (or to ulp(x) or ufp(x)-the definition of these functions is recalled at the end of this section). We will call these splittings relative splittings. In some languages, these splittings could be done using bit manipulations of the machine representations of the floating-point operands, but this would result in less portable programs.

The main contributions of this paper can be summarized as follows: we propose a new algorithm for computing floor(x) with rounding to nearest mode (Algorithm 2), a new FMAbased relative splitting with a smaller number of flops and latency than Veltkamp's splitting (Algorithm 4), and a survey of different algorithms to compute sign(x) • ulp(x) or sign(x) • ufp(x) up to a factor 2, some of which are underflow-safe and almost overflow-safe.

Interestingly, most of the algorithms presented in this paper have in common a sequence of two floating-point operations of the form "add some number z to the input, and then subtract z from the obtained sum."

In the following, we assume an IEEE-754 compliant floating-point arithmetic with radix β, precision p, and extremal exponents e min and e max . (In some cases we will only consider the case β = 2.) We denote by F the set of floatingpoint numbers in such an arithmetic. An element x of F can be written

x = M x • β ex-p+1 ,
where M x and e x are two integers such that |M x | < β p and e min e x e max , with |M x | maximum under these constraints. The significand of x is the number

M x • β -p+1 . A nonzero element x of F is normal if |x| β emin , and
subnormal otherwise. The absolute value of the significand of a normal floating-point number is at least 1 and less than β.

We will write RN to denote rounding to nearest for some given tie-breaking rule (assumed to be either "to even" or "to away", as in IEEE 754-2008). Hence, for t ∈ R and in the absence of overflow, RN(t) satisfies |t -RN(t)| |t -x| for all x ∈ F. We also write succ(t) to denote the successor in F of t (with the implicit assumption that if t is at least the largest element of F, then its successor should be +∞).

Before presenting splitting algorithms, let us recall some definitions of the functions ulp and ufp.

Definition I.1 (classical ulp). The unit in the last place of t ∈ R is

ulp(t) = β log β |t| -p+1 if |t| β emin , β emin-p+1 otherwise.
There are in fact some alternative definitions of ulp, which differ only around the powers of the radix and satisfy slightly different properties. In particular, Harrison's definition [START_REF] Harrison | A machine-checked theory of floating point arithmetic[END_REF] is as follows (and, to avoid confusion, we write ulp H ). Definition I.2 (Harrison's ulp). Let t ∈ R. The number ulp H (t) is the distance between the two closest straddling floating-point numbers x and y, so that x t y and x = y, assuming an unbounded exponent range.

One easily checks that in radix-2, precision-p arithmetic, if |t| is not a power of 2, then ulp(t) = ulp H (t), and if |t| = 2 k , then ulp(t) = 2 k-p+1 = 2ulp H (t), except in the subnormal range where ulp(t) = ulp H (t).

A related notion is that of ufp (unit in the first place), introduced by Rump, Ogita, and Oishi in [START_REF] Rump | Accurate floating-point summation, Part I: Faithful rounding[END_REF] and defined as follows.

Definition I.3 (ufp). The unit in the first place of t ∈ R is

ufp(t) = β log β |t| if t = 0, 0 if t = 0.

II. ABSOLUTE SPLITTINGS

The first algorithm we mention splits a floating-point number x into x h , which is x rounded to a nearest integer, and x = x -x h : The first occurrence we could find of this algorithm is in [5, p. 22], in radix 2 and with a constant C = 2 p-1 or C = 2 p-1 + 2 p-2 ; the latter constant aims at accommodating negative values of x and its use is referred to as "the 1.5 trick." A similar example is the code optimization guide of the AMD Athlon processor [6, p. 43], where the so-called "magic number" 6755399441055744, equal to 2 52 + 2 51 , is used to convert double-precision numbers to integers. It is also noted in [7, p. 5] that taking C = 2 52 yields a nearest integer function for double-precision numbers x ∈ [0, 2 48 ). In [START_REF] Rump | Accurate floating-point summation, Part I: Faithful rounding[END_REF], this algorithm is introduced (under the name ExtractScalar), still in radix 2 arithmetic, with the constant C = 2 k . In particular, it is shown in [START_REF] Rump | Accurate floating-point summation, Part I: Faithful rounding[END_REF] that if |x| 2 k , then x h is a multiple of 2 k-p (but needs not be a nearest integer when x is negative) and |x | 2 k-p .

Require: C, x ∈ F s ← RN(C + x) x h ← RN(s -C) x ← RN(x -x h ) {optional} return x h {or (x h , x )}
Of course, Algorithm 1 is nothing more than the so-called Fast2Sum algorithm applied to C and x, but with C considered as a given constant. As noted for example by Linnainmaa [START_REF] Linnainmaa | Analysis of some known methods of improving the accuracy of floating-point sums[END_REF], the sequence of operations of Fast2Sum already appears in early papers by Kahan [START_REF] Kahan | Further remarks on reducing truncation errors[END_REF] (for compensated summation) and Møller [START_REF] Møller | Quasi double-precision in floating point addition[END_REF] (with truncation instead of rounding to nearest) and was analyzed by Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], who proved that x is exactly the error of the addition of C and x when |C| |x| and β 3.

Theorem II.1. Assume p-digit floating-point numbers in radix β, and C integer with

β p-1 C β p . If β p-1 -C x β p -C, (1)
then the floating-point number x h returned by Algorithm 1 is an integer such that |x -

x h | 1/2 (that is, x h is equal to x rounded to a nearest integer). Furthermore, x = x h + x .
Proof. Adding C to Eq. ( 1) yields

β p-1 C + x β p , thus s = RN(C + x) lies in [β p-1 , β p ],
which implies s is an integer and |s -

(C + x)| 1/2. Now since C also lies in [β p-1 , β p ], s -C is exactly representable in precision p and thus x h = s -C is an integer: |x h -x| = |(s -C) -x| 1/2.
From x h = s -C and s = RN(C + x), we also deduce that x -x h equals C + x -RN(C + x), that is, the error of a floating-point addition. Hence x -x h must be a floating-point number (see for example [11, p. 236]) and we conclude that x = x -x h .

Note that Theorem II.1 scales in the obvious way: if

β k+p-1 C β k+p and β k+p-1 -C x β k+p -C, then x h is a multiple of β k , and |x | = |x -x h | 1 2 β k .
It is this property that is used in [START_REF] Rump | Accurate floating-point summation, Part I: Faithful rounding[END_REF] to "extract" parts of floatingpoint numbers that can be accumulated without error.

When x is nonnegative, the largest interval for x is attained for C minimum, i.e., C = β p-1 , and then the algorithm works for x β p -β p-1 . When x can be of either sign, the largest centered interval for x is attained when β p-1 -C and β p -C are opposite numbers, that is, C = (β p + β p-1 )/2, assuming β is even; then the algorithm works for |x| (β p-1 -β p )/2.

In particular, for β = 2 and p = 53, this yields the "magic number" C = 2 52 + 2 51 mentioned before, together with the range |x| 2 51 . Also, for β = 10 and p = 16 (decimal64 IEEE-754 format), the largest range for nonnegative inputs is x 9•10 15 , obtained with C = 10 15 ; and the largest range for signed inputs is |x| 4.5 • 10 15 , obtained with C = 5.5 • 10 15 .

A. Computing floor(x)

An interesting question is to compute x , or more generally x/2 k , for x a floating-point number. If the rounding mode is to nearest, how can we compute x efficiently?

Require: x ∈ F y ← RN(x -0.5) C ← RN(β p -x) s ← RN(C + y) x h ← RN(s -C) return x h Algorithm 2: Computing floor(x). Theorem II.2. Assume β is even, x ∈ F, 0 x β p-1 . Then Algorithm 2 returns x h = x .
Proof. We provide a detailed proof only for "ties to even", since similar arguments apply in the case where RN breaks ties "to away". Since β is even, 0.5 is in F. If 0 x < 0.5, then -0.5 y 0, C = β p , thus s = β p , and x h = 0 = x . Now assume 0.5 x β p-1 . In this case, it can be shown that x -0.5 is in F, so that

y = x -1/2. Write x = k + ε with k integer and 0 ε < 1; C is β p -k if ε < 1/2 or ε = 1/2
and k even, and β p -(k + 1) otherwise. In the case ε < 1/2 or ε = 1/2 and k even, we have s = RN(

β p + ε -0.5) = β p , thus x h = k. In the case ε > 1/2 or ε = 1/2 and k odd, we have s = RN(β p + ε -1.5) = β p -1, thus x h = RN((β p - 1) -(β p -(k + 1))) = k.

III. RELATIVE SPLITTINGS

We are now interested in expressing a precision-p floatingpoint number x as the exact sum of a (p -s)-digit number x h and an s-digit number x . The first use of such a splitting in the literature seems to be with s = p/2 : Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] used it to express the exact product of two floating-point numbers as a double-word. Another use is with s = p -1, so that x h is a power of β giving the order of magnitude of x. This can be desirable for at least two reasons:

• One may want to evaluate ulp(x) or ufp(x). These functions are very useful in the error analysis of floatingpoint algorithms. In that case, what we want to compute is fully determined (for example, for x nonzero, ufp(x) is the greatest integer power of β less than or equal to |x|). • One may just want to obtain an integer power of β close to |x|. Indeed, for scaling x, such a weaker condition will in general be enough, and one may hope that it can be satisfied using fewer operations than those needed to compute ufp(x).

A. Veltkamp's splitting and a new variant

We use here the description of Veltkamp's splitting from [START_REF] Linnainmaa | Software for doubled-precision floating-point computations[END_REF]. Given a floating-point number x of precision p in radix β, and an integer s < p, the following algorithm splits x into two non-overlapping floating-point numbers x h and x such that x = x h + x , with the significand of x h fitting in p -s digits, and the one of x in s digits (or even s -1 when β = 2 and s 2).

Require:

C = β s + 1 and x in F γ ← RN(Cx) δ ← RN(x -γ) x h ← RN(γ + δ) x ← RN(x -x h ) return (x h , x )
Algorithm 3: Veltkamp's splitting.

Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] analyzes this algorithm in radix 2, with the implicit assumption that no overflow occurs. This analysis was extended to any radix β by Linnainmaa in [START_REF] Linnainmaa | Software for doubled-precision floating-point computations[END_REF]. Note that the algorithm works correctly even in the presence of underflows (due to the special shape of C and since underflowing additions are exact, as noted for example in [START_REF] Hauser | Handling floating-point exceptions in numeric programs[END_REF]). Moreover, Boldo [START_REF] Boldo | Pitfalls of a full floating-point proof: Example on the formal proof of the Veltkamp/Dekker algorithms[END_REF] shows for any radix β that if Cx does not overflow, then no other operation will overflow.

If an FMA instruction is available, we suggest to use the following variant, that requires fewer operations.

Require:

C = β s + 1 and x in F γ ← RN(Cx) x h ← RN(γ -β s x) x ← RN(x -x h ) {or x ← RN(Cx -γ)} return (x h , x )
Algorithm 4: FMA-based relative splitting.

Theorem III.1. Let x ∈ F and s ∈ Z such that 1 s < p.

Then, barring underflow and overflow, Algorithm 4 computes x h , x ∈ F such that x = x h + x and, if β = 2, the significands of x h and x have at most p -s and s bits, respectively. If β > 2 then their significands have at most p -s + 1 and s + 1 digits, respectively.

Proof. The result is clear for x = 0 and, for x nonzero, we can restrict with no loss of generality to 1 x < β. For s 0, this implies β s Cx < β s+1 + β β s+2 , so ulp(Cx) β s-p+1 and thus γ = RN(Cx) is a multiple of β s-p+1 . Since x 1, β s x is also a multiple of β s-p+1 , and thus

γ -β s x = M • β s-p+1
for some integer M . Now, let be such that γ = Cx + . We have | |

1 2 β s-p+2 since ulp(Cx) β s-p+2
, and, on the other hand, γ -β s x = x + .

Therefore, recalling that |x| < β,

|M | |x| + | | β s-p+1 < β p-s + 1 2 β. (2)
This implies |M | β p and that γ -β s x is in F. Consequently, x h = x + and x = RN(-). Since x h is a multiple of β s-p+1 and x is a multiple of β 1-p , their difference is also a multiple of β 1-p and, using the bound on | | seen above, there exists an integer E such that

= E • β 1-p , |E| 1 2 β s+1 . (3)
Hence |E| β p and thus ∈ F and x = x -x h , as wanted. From ( 2) and (3) we deduce that |M | β p-s+1 and |E| β s+1 , which means that the significands of x h and x fit into p -s + 1 and s + 1 digits, respectively. If β = 2 then the strict inequality in (2) implies |M | 2 p-s , and (3) gives |E| 2 s . Hence p -s and s bits suffice to represent the significands of x h and x , respectively, in this case.

Finally, note that since x h = x + and γ = Cx + , we have

x -x h = Cx -γ, so the suggested variant for x can be used as well.

We conclude this section with several remarks about Algorithm 4:

• If an efficient FMA is available, it can be used for the computation of x h , and Algorithm 4 takes only 3 flops instead of 4 for Algorithm 3. In fact, if for x we use the variant x ← RN(Cx -γ), then x can be obtained simultaneously with x h , and this reduces the depth (or latency) of the relative splitting from 4 to 2. • Even if an efficient FMA is not available, the two floatingpoint numbers γ and β s x can be computed in parallel, and so, again, the depth of Veltkamp's splitting is reduced, but only from 4 to 3. • If β > 2, the bounds p -s + 1 and s + 1 on the numbers of digits needed to represent the significands of x h and x can be attained. For example, if x = β -β 1-p , then it is easy to check that γ = β s+1 + β and

x h = (1 + β s-p ) • β, x = -(β s + 1) • β 1-p .
• If β = 2, the bounds p -s and s specific to this case are attainable as well. To show this, we can consider the same example as before, namely x = 2-2 1-p , and check that we then have γ = (1 + 2 -s -2 1-p ) • 2 s+1 and

x h = (2 p-s -1) • 2 s-p+1 , x = (2 s -1) • 2 1-p .

B. A special case: extracting one bit only

Assuming radix 2, let us now consider the problem of computing ufp(x) or ulp(x), and more generally the question of scaling x, that is, deducing from x an integer power of 2 that is "close" to |x|.

For scaling, a possible solution is to determine a nearest power of two. This can be done by applying Veltkamp's splitting (Algorithm 3) to x with s = p -1: in this case, the resulting x h has a 1-bit significand and we know from [14, Theorem 1] that it is nearest x in precision p -s = 1.

For computing sign(x) • ufp(x), we can use the following algorithm, introduced by Rump in [START_REF] Rump | Ultimately fast accurate summation[END_REF].

Require: β = 2, ϕ = 2 p-1 + 1, ψ = 1 -2 -p , and x ∈ F q ← RN(ϕx) r ← RN(ψq) δ ← RN(q -r) return δ
Algorithm 5: Computing sign(x)•ufp(x) for radix 2, as in [15, Algorithm 3.5].

These solutions, however, raise the following issues.

• If |x| is large, then an overflow can occur in the first line of both Algorithms 3 and 5. It is always possible to slightly change Algorithm 5, by scaling it, in order to avoid overflow: it suffices for instance to replace ϕ by 1 2 +2 -p and to return 2 p δ at the end. However, this variant will not work for subnormal x. Thus, to use Algorithm 5, we somehow need to check the order of magnitude of x.

β k x β k+1 -β k-p+1 : 1) if x = β k then δ = x; 2) if β k + β k-p+1 x < β 2 • β k then ulp(ψ * x) = β k+1 |ψ * x -β p x| = |x| < 1 2 β k+1 , so that RN(ψ * x) = β p x and δ = 0; 3) if x > β 2 • β k , we still have ulp(ψ * x) = β k+1
, and we deduce that

β p x -ulp(β p x) ψ * x = β p x -x < β p x - 1 2 ulp(β p x);
therefore, RN(ψ * x) = β p x -β k+1 , so that δ = β k+1 > x; 4) the case x = β 2 • β k (which can occur only when the radix β is even) is similar to case 2) or case 3), depending on the way RN breaks ties.

C. Alternative algorithms for computing sign(x) • ulp(x) and sign(x) • ufp(x) If we replace the constant ϕ in Algorithm 5 by ϕ/2 p-1 and if we ignore the possibility of underflow and overflow, we obtain a way to compute sign(x) • ulp(x) which uses three floating-point operations.

An alternative method is given by Algorithm 8 below, which, if an FMA is available, uses only two floating-point operations. This method uses the same scheme as [16, Algorithm 2], outside the subnormal range, except that in [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF] the constant ψ = u(1 + 2u) = succ(u) is used, whereas we use here 3 2 u, with u = 2 -p . We determined experimentally that, with or without FMA, ψ = u(1 + 2u) = succ(u) is the smallest constant that works, and ψ = succ( 3 2 u) is the largest one.

Require: β = 2, ψ = 2 -p + 2 -p-1 , and x ∈ F \ {0} a ← RN(x + ψx) {or RN(x + RN(ψx)) without FMA} δ ← RN(a -x) return δ
Algorithm 8: Computing sign(x)•ulp(x) for radix 2 and x = 0.

Theorem III.3. For x = 0 and p 2, if no underflow or overflow occurs, then Algorithm 8 returns sign(x) • ulp(x).

Proof. For x < 0, it is clear that the algorithm returns the opposite of what it returns for -x. Let us now assume x > 0.

Without loss of generality we can assume 1 x < 2. With u = 1 2 ulp(1) = 2 -p , note that ψ = 3 2 u. For x = 1, ψx = ψ is exactly representable, thus with or without FMA, we get a = RN(1 + ψ) = 1 + 2u, since the rounding boundaries for 1 + 2u are 1 + u and 1 + 3u.

It remains to prove the theorem for 1 + 2u x 2 -2u. We then have

3 2 u(1 + 2u) ψx 3 2 u(2 -2u).
It follows that x + u < x + ψx < x + 3u, thus with an FMA, x + ψx is rounded to x + 2u, and thus δ = 2u. Now consider the case without an FMA. On the left-hand side, since 3 2 u is exactly representable, because p 2, we have 3 2 u RN(ψx). On the right-hand side, 3 2 u(2 -2u) is strictly less than the midpoint u(3 -2u), so that RN(ψx) u(3 -4u) < 3u. Since x + u < x + RN(ψx) < x + 3u, x + RN(ψx) is rounded to x + 2u, and δ = 2u.

To compute sign(x)•ufp(x), one can simply scale x by 2 p-1 within Algorithm 8, which can be performed with 3 FMA instead of 2, but still with a depth of 2: we first compute a = RN(x + ψx) and y = RN(2 p-1 x) in parallel, and then RN(2 p-1 a -y) yields sign(x) • ufp(x).

It does not seem possible to compute sign(x) • ufp(x) with only 2 FMA. However, perhaps surprisingly, 2 FMA suffice to obtain sign(x) • 4ufp(x):

a ← RN((2 p+1 -2) • x -x) δ ← RN(2 p+1 • x -a)
Scaling this last operation by 2 -p-1 yields the following variant of Algorithm 8 for computing sign(x) • ulp(x), still with 2 FMA:

a ← RN((2 p+1 -2) • x -x) δ ← RN(x -2 -p-1 • a)
IV. UNDERFLOW-SAFE AND ALMOST OVERFLOW-FREE

SCALING

In this section we assume that β = 2 and p 4 (and we recall that RN breaks ties "to even" or "to away", as defined in IEEE 754-2008). Furthermore, we denote by η the smallest positive element of F:

η = 2 emin-p+1 .
Given a nonzero floating-point number x, we want to compute a scaling factor δ that satisfies the following properties:

• |x|/δ is much above the underflow threshold and much below the overflow threshold (so that, for example, we can safely square it); • δ is an integer power of 2 (in order to avoid rounding errors when multiplying or dividing by it). The algorithms proposed in Section III-B are simple, but underflow or overflow can occur for many inputs x, leading to an output δ that does not satisfy the above properties. Let us now introduce an algorithm that is underflow-safe and almost overflow-free in the sense that only the two extreme values x = ±(2 -2 1-p ) • 2 emax must be excluded.

The first three lines of Algorithm 9 mimic Algorithm 1 of [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF], with slightly different notation.

A. An example of application

Let us consider the use of Algorithm 9 for computing √ a 2 + b 2 . Assume we wish to evaluate √ a 2 + b 2 without risking spurious overflows and underflows. Assume RN is with "ties to even" (so that the scaling factor provided by Algorithm 9 is, up to two extreme values, always a power of 2) and, in addition, that e min < -p,

(25/2) • 2 2p 2 emax+1 , (4)
which always holds with the usual binary floating-point formats.

In the following, we assume that once the scaling is performed, we just use the naive formula. (Although other evaluation schemes are possible, we do not deal with them in this paper, which focuses only on the scaling part).

Let

c = RN RN a 2 + RN b 2 .
The computation of c is obviously overflow-free. We will deal with possible loss of precision due to underflow later on. Define δ c as the value returned by Algorithm 9 with c as input. We suggest to use δ c as scaling factor. As a consequence, since we use the naive formula with the scaled inputs, we actually compute

Since δ c is a power of 2-we assumed RN is with "ties to even", in which case the situation where δ c = 3ulp(c) cannot happen-, multiplication and division by δ c can induce a rounding error only when the result is a subnormal number. Let us quickly examine this scaling choice. Without loss of generality, we can assume that a 0 and a |b|. Note that if a 2 emin then RN(a/2) = a/2, and that if 2η a 2 emin , we have |RN(a/2) -a/2| η/2. From this we deduce that if 2η a we always have

2 5 a RN a 2 2 3 a. (6)
(The lower bound is attained for a = 5η, and the upper bound it attained for a = 3η.) 1) If RN(a/2) = 0 (which implies a 2η) then RN(a/2) c 2RN(a/2), so that, using (6), 2a/5 c 4a/3. Therefore, from Theorem IV.2, 3 4

3 4 • c δ c a δ c 5 2 • c δ c 5 2 • (2 p -1).
Hence, no underflow can result from the computation of RN(a/δ c ) 2 and, using (4), no overflow can result from the computation of RN(RN(a/δ c ) 2 + RN(b/δ c ) 2 ).

If an underflow results from the computation of (b/δ c ) 2 (i.e., if |b|/δ c < 2 emin/2 ) then one can show using (4) that RN(RN(a/δ c ) 2 + RN(b/δ c ) 2 ) = RN(a/δ c ) 2 , which implies that the calculation in (5) will return a (see for example [START_REF] Kahan | Mathematics written in sand-the HP-15C, Intel 8087, etc[END_REF]), and that this is a very accurate result. If no underflow occurs then the scaling induces no error. As a consequence, if RN(a/2) = 0 the usual relative error bound 2 -p+1 holds unless the final result is subnormal.1 2) If RN(a/2) = 0 then a and b are 0 or ±η, and c = 0. In such a case, we take δ c = η, which is the best possible splitting constant: the final result will be 0 if a = b = 0, and η otherwise, i.e., the correct rounding of √ a 2 + b 2 .

V. EXPERIMENTAL RESULTS

In this section we compare our algorithms with bit manipulation algorithms, and with standard functions of the mathematical library, in the C language. Our experimental framework consists of an Intel i5-4590 processor running at 3.3GHz, under Debian testing (buster), with gcc 7.3.0, and optimization level -O3. We also set the FPU control word so that all roundings are done in double precision (and not in double-extended precision as by default). We first generate an array of 10 8 random binary64 numbers in the range [0, 2 32 ], and we then compute the time to call each function on all numbers. The C programs we used are available at http://homepages.loria.fr/PZimmermann/papers/#split. round floor Algorithms 1 and 2 0.106s 0.173s Bit manipulation 0.302s 0.203s GNU libm 0.146s 0.209s Fig. 1. Rounding to an integer. Note: we used the GNU libm rint function for rounding to nearest (which is faster than nearbyint), and floor for rounding towards -∞.

x h x time Algorithm 3 26 bits 26 bits 0.108s Algorithm 4 26 bits 27 bits 0.106s Algorithm 4 with FMA 26 bits 27 bits 0.108s Bit manipulation 26 bits 27 bits 0.095s These results show that the algorithms we presented -in their respective domain of validity -are competitive with the less portable bit manipulation algorithms, and in some cases faster. Even, in some languages, for example in Javascript, there is no (efficient) way to manipulate the bits of a floatingpoint number, thus the user is left with no alternative.

VI. CONCLUSION AND FUTURE WORK

We have presented several splitting algorithms that use floating-point arithmetic with rounding to nearest. Some have been part of the floating-point folklore for decades, and some are new. These algorithms have many applications: we show that, beyond the traditional contexts of extended-precision multiplication, accurate summation, and computing the ulp and ufp functions, floating-point splittings can also be used to implement some integer functions such as floor or round, as well as some scalings aimed at reducing the risk of spurious underflows and overflows. In radix-2 floating-point arithmetic, simple algorithms can evaluate sign(x) • ulp(x), sign(x) • ulp H (x), and sign(x) • ufp(x), and also compute convenient scaling factors. However, generalizing this property to radix β (or, at least, to radix 10) does not seem to be straightforward and will have to be investigated in a future work.

Algorithm 1 :

 1 Rounding to a nearest integer, under conditions specified by Theorem II.1. The input value is x, and the algorithm uses a constant C.

  RN δ c × RN RN (RN (a 2 ) + RN (b 2 )) , with a = RN a δc and b = RN b δc .

Fig. 2 .

 2 Fig. 2. Splitting a double-precision number into x h and x .

Fig. 3 .

 3 Fig. 3. Comparison of different ways of computing sign(x) • ufp(x) (possibly up to a factor 2 for Algorithm 7).

•

  If we are only interested in scaling x, then requiring to produce the exact value of ufp(x) may seem like overkill: maybe one can get sometimes ufp(x) and sometimes Based on these remarks, we suggest considering the following algorithms as well. (Recall that ulp H refers to Harrison's ulp as in Definition I.2.) Algorithm 6 uses only 2 floating-point operations, and if an FMA is available, this also applies to Algorithm 7. In the absence of an FMA, 3 operations are needed for Algorithm 7, similarly to Algorithm 5. • Both Algorithms 4 and 7 are based on a sequence of operations of the form RN RN (C + 1)x -RN(Cx) for some prescribed constant C. Specifically, C = β s in Algorithm 4 and C = 2 p -1 in Algorithm 7. Note also that Algorithm 3.6 in [2] is quite similar to Algorithm 7, except that it uses the constant C = 2 p instead of C = 2 p -1, which requires a special treatment when x is a power of 2. • Note finally that by Theorem III.2, the output of Algorithm 7 is equal to x if and only if |x| is a power of two. Algorithm 7 does not straightforwardly generalize to higher radices: if one sets ψ * to β p -1 and if δ becomes RN(β p x -a), then we do not always have δ equal to a power of β close to x. However, we still have the property that δ = x if and only if x is a power of the radix. Hence Algorithm 7 provides a simple way to check this property using only floating-point operations. More precisely, assuming

	1 2 ufp(x) or 2ufp(x) (or any other convenient value) with
	a cheaper algorithm.
	Theorem III.2. If |x| > 2 emin , then Algorithm 6 returns
	sign(x)•ulp H (x) = sign(x)•	1 2 ulp(x) if |x| is a power of 2, ulp(x) otherwise.
	If |x| < 2 emax-p+1 , then Algorithm 7 returns
	sign(x)•	ufp(x) if |x| is a power of 2, 2ufp(x) otherwise.

Require: β = 2, ψ = 1 -2 -p , and x ∈ F a ← RN(ψx) δ ← RN(x -a) return δ Algorithm 6: Computing sign(x) • ulp H (x) for radix 2 and |x| "large enough" (namely, |x| > 2 emin ).

Require:

β = 2, ψ * = 2 p -1, and x ∈ F a ← RN(ψ * x) δ ← RN(2 p x -a) return δ

Algorithm 7: Computing sign(x) • ufp(x) up to a factor 2, for radix 2 and |x| "small enough" (namely, |x| < 2 emax-p+1 ).

Proof. The first claim is a direct consequence of

[START_REF] Rump | Ultimately fast accurate summation[END_REF] Lemma 3.6

]. The second claim then follows from the fact that Algorithm 7 is Algorithm 6 applied to 2 p x.

Again, let us conclude with several comments.

•

If the final result is subnormal (i.e., if an underflow happens during the final multiplication by δ: it cannot happen before), then we need to add an additional absolute error η/2 due to the final rounding, so that we end up with an absolute error bound 3η/2.

ACKNOWLEDGMENT

We thank Pierrick Gaudry for the initial question that led to that work, and for comments on a first version of this paper. This work was partially supported by the French National Research Agency grants ANR-13-INSE-0007 (MetaLibm project) and ANR-14-CE25-0018-01 (FastRelax project).

Require: β = 2, Φ = 2 -p + 2 -2p+1 , η = 2 emin-p+1 , and x ∈ F. y ← |x| e ← RN(Φy + η)

{or e ← RN(RN(Φy) + η) without FMA} y sup ← RN(y + e) δ ← RN(y sup -y) return δ Algorithm 9: Computing a scaling factor that is overflow-free and underflow-safe, for radix 2 and |x| = (2 -2 1-p ) • 2 emax .

Below we recall a subcase of a result obtained by Rump et al. [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF].

Theorem 2.2 in [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF] was proven for e equal to RN(RN(Φ|x|) + η) (that is, for the "without FMA" version of Algorithm 9). However, one can easily check that the proof given in [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF] still holds for e = RN(Φ|x| + η) (the "with FMA" version). An immediate consequence of Theorem IV.

Theorem 2.2 from [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF] has been formally checked using the Coq automatic proof checker, which gives much confidence in this result. However, in [START_REF] Rump | Computing predecessor and successor in rounding to nearest[END_REF] it is claimed further that if |x| ∈ [2 emin , 2 emin+2 ] then y sup is succ(succ(|x|)). As we are going to see, this property does not always hold. Let y = |x| and consider the following cases:

1) If 2 emin y 2 emin+1 -2η then ulp(x) = η. A simple calculation shows that

so that e = RN(Φy + η) = 2η = 2ulp(x). An immediate consequence is that RN(y + e) = y + e = succ(succ(|x|)). Given the domain of x, we deduce that

, so e = RN(Φy + η) = 2η = 2ulp(x).

Therefore y + e is equal to the midpoint 2 emin+1 + η and we have the following:

• if RN is with "ties to even" then RN(y + e) = 2 emin+1 = succ(y), so δ = η = ulp(x);

• if RN is with "ties to away" then RN(y + e) = 2 emin+1 + 2η = succ(succ(y)), so δ = 3η = 3ulp(x).

In this case, δ is not a power of 2. As we will see this is the only case where this happens; 3) If 2 emin+1 y 2 emin+2 -2η then ulp(x) = 2η. A simple calculation shows that 2η+2 emin-2p+2 Φy+η 3η+2 emin-2p+2 -2 emin-3p+3 so that e = RN(Φy + η) ∈ {2η, 3η} (and one easily checks that both cases actually do occur). When e = 2η, we obtain y sup = succ(y) and

When e = 3η, depending on the tie-breaking rule and (for ties-to-even) on whether the last bit of y is a zero or a one, we obtain either y sup = succ(y) and Consequently, depending on the tie-breaking rule, if |x| ∈ [2 emin , 2 emin+2 ] then Algorithm 9 will return ulp(x) or 2ulp(x) or, in one case only, 3ulp(x). We deduce the following theorem.

Theorem IV.2. For x ∈ F with |x| = (2 -2 1-p ) • 2 emax , the value δ returned by Algorithm 9 satisfies the following:

• if RN is with "ties to even" then δ is a power of 2;

• if RN is with "ties to away" then δ is a power of 2, unless |x| = 2 emin+1 -2 emin-p+1 , in which case it equals

Therefore, Algorithm 9 is a very good candidate for performing scalings.

Interestingly enough, in single precision/binary32 arithmetic (e min = -126, p = 24), the number (2 emin+1 -2 emin-p+1 )/(3• 2 emin-p+1 ) is a floating-point number, so that even for RN with "ties to away" and |x| = 2 emin+1 -2 emin-p+1 , the calculation of x/δ is error free. This is not the case in double precision/binary64 arithmetic.


Xi Wu
email: xiwu@sei.ecnu.edu.cn

Ximeng Li

Alberto Lluch Lafuente

Flemming Nielson

Hanne Riis Nielson

Klaim-DB: A Modeling Language for Distributed Database Applications

We present the modelling language, Klaim-DB, for distributed database applications. Klaim-DB borrows the distributed nets of the coordination language Klaim but essentially re-incarnates the tuple spaces of Klaim as databases, and provides high-level language abstractions for the access and manipulation of structured data, with integrity and atomicity considerations. We present the formal semantics of Klaim-DB and illustrate the use of the language in a scenario where the sales from different branches of a chain of department stores are aggregated from their local databases. It can be seen that raising the abstraction level and encapsulating integrity checks (concerning the schema of tables, etc.) in the language primitives for database operations benefit the modelling task considerably.

Introduction

Today's data-intensive applications are becoming increasingly distributed. Multinational collaborations on science, economics, military etc., require the communication and aggregation of data extracted from databases that are geographically dispersed. Distributed applications including websites frequently adopt the Model-View-Controller (MVC) design pattern in which the "model" layer is a database. Fault tolerance and recovery in databases also favors the employment of distribution and replication. The programmers of distributed database applications are faced with not only the challenge of writing good queries, but also that of dealing with the coordination of widely distributed databases. It is commonly accepted in the formal methods community that the modelling of complex systems in design can reduce implementation errors considerably [START_REF] Abrial | Formal methods: Theory becoming practice[END_REF][START_REF] Woodcock | Formal methods: Practice and experience[END_REF].

Klaim [START_REF] De | KLAIM: A kernel language for agents interaction and mobility[END_REF] is a kernel language for specifying distributed and coordinated processes. In Klaim, processes and information repositories exist at different localities. The information repositories are tuple spaces, that can hold data and code. Processes can read tuples from (resp. write tuples to) local or remote Part of the work was done when Xi Wu was a visiting researcher at DTU compute.

tuple spaces, or spawn other processes to be executed at certain localities. Many distributed programming paradigms can be modelled in Klaim.

While Klaim provides an ideal ground for the modelling of networked applications in general, the unstructured tuple spaces and low-level operations mostly targeting individual tuples create difficulty in the description of the datamanipulation tasks usually performed using a high-level language such as SQL. A considerable amount of meta-data needed by databases has to be maintained as particular tuples or components of tuples, the sanity checks associated with database operations have to be borne in mind and performed manually by the programmer, difficulties arise when batch operations are performed and atomicity guarantees are needed, and so on.

To support the modelling of applications operating on distributed, structured data, we propose the language Klaim-DB, which is inspired by Klaim in that it allows the distribution of nodes, and remote operations on data. Its succinct syntax eases the precise formulation of an operational semantics, giving rigor to high-level formal specification and reasoning of distributed database applications. The language also borrows from SQL, in that it provides structured data organized as databases and tables, and high-level actions that accomplish the data-definition and data-maniputation tasks ubiquitous in these applications. We use a running example of database operations in the management of a large-scale chain of department stores. Each individual store in the chain has its local database containing information about the current stock and sales of each kind of product. The semantic rules for the core database operations will be illustrated by the local query and maintenance of these individual databases, and our final case study will be concerned with data aggregation across multiple local databases needed to generate statistics on the overall product sales.

This paper is structured as follows. In Section 2, the syntax of Klaim-DB is presented, which is followed by the structural operational semantics specified in Section 3, with illustration of the semantic rules for the main database operations. Our case study is then presented in Section 4. Extensions of Klaim-DB and a discussion of alternative language design choices are covered in Section 5. We conclude in Section 6 with a discussion of related works.

Syntax

The syntax of Klaim-DB is presented in Table 1. A net N models a database system that may contain several databases situated at different localities. As in standard Klaim [START_REF] De | KLAIM: A kernel language for agents interaction and mobility[END_REF], we distinguish between physical localities, also called sites (s ∈ S), and logical localities ( ∈ L) that are symbolic names used to reference sites. At each site s, there is an "allocation environment" ρ : L → S mapping the logical localities known at s to the sites referred to by them.

We assume for simplicity that there is only one database at each site. The syntax s :: ρ C for a node of the net captures the ensemble C of processes and tables of the database at the site s, where the allocation environment is ρ. The parallel composition of different nodes is represented using the || operator. With

::= nil | N1||N2 | (νs)N | s ::ρ C C ::= P | (I, R) | C1|C2 P ::= nil | a.P | A(ẽ) | foreachs T in R : P | foreachp T in R : P | P1; P2 a ::= insert(t, tb)@ | insert tb(TBV , tb)@ | delete(T, ψ, tb, !TBV )@ | sel ext(T, ψ, tb, t, !TBV )@ | sel int(T, ψ, TBV , t, !TBV ) | update(T, ψ, t, tb)@ | aggr(T, ψ, tb, f, T )@ | create(I)@ | drop(tb)@ t ::= e | t1, t2 T ::= e | !x | T1, T2
the restriction operator (νs)N , the scope of s is restricted to N . Each table is represented by a pair (I, R) where I is an interface that publishes certain structural information about the table as attributes (for example, I.id stands for the table identifier and I.sk is a schema describing the data format of the table), and R is a multiset of tuples, in which each tuple corresponds to one row of data. The construct C 1 |C 2 is the parallel composition of different components.

We distinguish between tuples t and templates T . A template T can contain not only actually fields that are expressions e, but also formal fields !x where x is a variable that can be bound to values.

A process P can be an inert process nil, an action prefixing a.P , a parameterized process invocation A(ẽ), a looping process foreach s T in R : P , or foreach p T in R : P , or a sequential composition P 1 ; P 2 , for which we require that bv (P 1 ) ∩ fv (P 2 ) = ∅. Looping is introduced in addition to recursion via process invocation, to ease the task of traversing tables or data selection results in a natural way. We also allow both prefixing and sequential composition in our language, as in CSP [START_REF] Hoare | Communicating Sequential Processes[END_REF]. Sequential composition is needed to facilitate the specification of possible continuation after the completion of a loop, whereas prefixing is retained and used in situations where substitutions need to be applied after receiving an input. The difference between the two variants of looping process is that the sequential loop foreach s T in R : P goes through the different rounds sequentially, while the parallel loop foreach p T in R : P forks one parallel process for each round.

A process can perform nine different kinds of actions. Actions insert(t, tb)@ , insert tb(TBV , tb)@ , delete(T, ψ, tb, !TBV )@ , sel ext(T, ψ, tb, t, !TBV )@ , sel int (T, ψ, TBV , t, !TBV ), update( T, ψ, t, tb)@ and aggr(T, ψ, tb, f, T )@ are used to access/manipulate the data inside a table; they resemble the operations performed by the "Data-Manipulation Language" in SQL. On the other hand, actions create(I)@ , and drop(tb)@ are used for the creation and deletion of a table -they correspond to the operations performed by the "Data-Definition Language" in SQL.

The actions insert(t, tb)@ and insert tb(TBV , tb)@ are used to insert a new row t, or all the rows of a table bound to the table variable TBV , into a table named tb inside the database at , respectively. On the other hand, the action delete(T, ψ, tb, !TBV )@ deletes all rows matching the pattern T and the predicate ψ from table tb in the database located at , and binds the deleted rows to the table variable TBV .

The language has two variants of "selection": an "external" one that selects data from tables actually existing in databases, identified by their table identifiers, and an "internal" one that selects data from temporary tables bound to table variables. The action sel ext(T, ψ, tb, t, !TBV )@ performs the "external" selection. It picks all rows matching the pattern T as well as satisfying the predicate ψ, from the table identified by tb of the database located at , and binds the resulting table into the table variable TBV . On the other hand, sel int(T, ψ, TBV , t, !TBV ) performs the "internal" selection, i.e., from the content of the table variable TBV , and binds the resulting table further into the table variable TBV . In sel int(T, ψ, TBV , t, !TBV ), the meanings of T , ψ, and t are the same as those in sel ext(T, ψ, tb, t, !TBV )@ . In both variants, we require that each component of t should be a value or a bound variable of T .

The action update(T, ψ, t, tb)@ replaces each row matching T yet satisfying ψ in table tb (at ) with a new row t, while leaving the rest of the rows unchanged. It is required that fv (t) ⊆ bv (T ).

The action aggr(T, ψ, tb, f, T )@ applies the aggregation function f on the multiset of all rows meeting T and ψ in table tb (at ) and binds the aggregation result to the pattern T .

The action create(I)@ (resp. drop(tb)@ ) creates a table identified by I (resp. drops the table identified by tb) in the database at . An item is a row in a table containing sequences of values, which can be obtained from the evaluation of an expression e. Pattern-matching is used to manage the data inside a table by means of a given template T , which is a sequence of values and variables.

Setting the Scene for the Running Example Consider the management of a large-scale chain of department stores, in which the head office can manage the sales of different imaginative brands (e.g., KLD, SH,...) in its individual branches, as shown in Figure 1.

We will use underlined symbols for logical and physical localities, as well as allocation environments, to distinguish between their uses in the example and elsewhere (e.g., in the semantics). Suppose the database of the head office is maintained on machine s 0 and the databases of its branches are maintained on machines s 1 to s n . The site s 0 has the local environment

ρ 0 = [self → s 0 ][ 1 → s 1 ]...[ n → s n ]. We use ρ j = [self → s j ]

Fig. 2. The Table Stores

A table with identifier Stores exists in the database of the head office, and records the information of its branches, as shown in Figure 2. The header partially describes the schema I 0 of the table, and the subsequent rows, constituting the multiset R 0 , contain the information of the different branches. Formally, we have I 0 .id = Stores, I 0 .sk =< "City" : String, "Address" : String, "Shop name" : String, "Brand" : Set, "Logical Locality" : String >.

Each database of a branch has several tables identified by the name of the brand, which record the information of the stock and sales of the corresponding brand. The table of one shoe brand, KLD, in one branch, is shown in Figure 3. The identifier of this table is KLD and the schema is < "Shoe ID" : String, "Shoe name" : String, "Y ear" : String, "Color" : String, "Size" : Int, "In-stock" : Int, "Sales" : Int >.

Fig. 3. The Table KLD in one branch

To sum up, the network of databases and operating processes can be represented by

s0 ::ρ 0 ((I0, R0)|C 0 ) || s1 ::ρ 1 ((I1, R1)|C 1 ) || ... || sn ::ρ n ((In, Rn)|C n ),
where for j ∈ {1, ..., n}, (I j , R j ) describes the local table for the brand KLD inside its database at s j , and for each k ∈ {0, ..., n}, C k stands for the remaining processes and tables at the site s k .

Semantics

We devise a structural operational semantics [START_REF] Plotkin | A structural approach to operational semantics[END_REF] for Klaim-DB, as shown in Table 3 andTable 4. The semantics is defined with the help of a structural congruence -the smallest congruence relation satisfying the rules in Table 2, where the α-equivalence of N and N is denoted by

N ≡ α N . Table 2. The Structural Congruence N1||N2 ≡ N2||N1 (νs1)(νs2)N ≡ (νs2)(νs1)N (N1||N2)||N3 ≡ N1||(N2||N3) N1||(νs)N2 ≡ (νs)N1||N2 (if s ∈ fn(N1)) N ||nil ≡ N s ::ρ C ≡ s ::ρ C|nil N ≡ N (if N ≡α N ) s ::ρ (C1|C2) ≡ s ::ρ C1 || s ::ρ C2
We start by explaining our notation used in the semantic rules. The evaluation function E • ρ evaluates tuples and templates under the allocation environment ρ. The evaluation is performed pointwise on all the components of tuples t and templates T , much in the manner of standard Klaim [START_REF] De | KLAIM: A kernel language for agents interaction and mobility[END_REF]. We denote by [a/b] the substitution of a for b; in particular, when a and b are tuples of the same length, point-wise substitution is represented. Pattern matching is performed with the help of the predicate match(eT, et) where eT is an evaluated template and et is an evaluated tuple. This match(eT, et) can be defined in a manner much like that of [START_REF] De | KLAIM: A kernel language for agents interaction and mobility[END_REF]. In more detail, match(eT, et) gives a boolean value indicating whether the pattern matching succeeded. In case it did, the substitution σ = [et/eT ] can be resulted. For example, match((3, !x, !y), (3, 5, 7)) = true and [(3, 5, 7)/(3, !x, !y)] = [5/x, 7/y]. We denote by ψσ the fact that the boolean formula ψ is satisfied after the substitution σ is applied to it. We will use , ∩ and \ to represent the union, intersection and substraction, of multisets and the detailed definitions are given in Appendix A.

The notation I[a → b][...] represents the update of the interface I that maps its attribute a to b; multiple updates are allowed. In addition, we will use I.sk ↓ T t to represent the projection of the schema I.sk according to the template T (matching the format requirements imposed by I.sk ) and the tuple t. When t only contains constants or the bound variables of T as its components, I.sk ↓ T t is a new schema that describes only the columns referred to in the variable components of t. Since we have left the schema under-specified, this projection operation is illustrated in Example 1, rather than formally defined.

Example 1. Suppose I.sk =< "Shoe ID" : String, "Shoe name" : String, "Y ear" : String, "Color" : String, "Size" : Int, "In-stock" : Int, "Sales" : Int >, which specifies that the table having I as its interface has seven columns: "Shoe ID of type String, "Shoe name" of type String, "Year" of type String, "Color" of type String, "Size" of type Int, "In-stock" of type Int, "Sales" of type Int. Suppose in addition that T = ("001", "HighBoot", !x, !y, !z, !w, !p) and t = (y, z, p). Then I.sk ↓ T t =< "Color" : String, "Size" : Int, "Sales" : Int >.

ρ1(l) = s2 I.id = tb t |= I.sk R = R {E t ρ 1 } s1 ::ρ 1 insert(t, tb)@ .P || s2 ::ρ 2 (I, R) → s1 ::ρ 1 P || s2 ::ρ 2 (I, R ) (INS TB) ρ1(l) = s2 I.id = tb I .sk = I.sk R = R R s1 ::ρ 1 insert tb((I , R ), tb)@ .P || s2 ::ρ 2 (I, R) → s1 ::ρ 1 P || s2 ::ρ 2 (I, R ) (DEL) R = {t | t ∈ R ∧ ¬(match(E T ρ 1 , t) ∧ ψ[t/E T ρ 1 ])} R = {t | t ∈ R ∧ match(E T ρ 1 , t) ∧ ψ[t/E T ρ 1 ]} ρ1(l) = s2 I.id = tb σ = [(I, R )/TBV ] s1 ::ρ 1 delete(T, ψ, tb, !TBV )@ .P || s2 ::ρ 2 (I, R) → s1 ::ρ 1 P σ || s2 ::ρ 2 (I, R ) (SEL EXT) ρ1(l) = s2 I.id = tb I = I[id → ⊥][sk → I.sk ↓ T t ] σ = [(I , R )/TBV ] R = {E tσ ρ 1 | ∃t : t ∈ R ∧ match(E T ρ 1 , t ) ∧ σ = [t /E T ρ 1 ] ∧ ψσ} s1 ::ρ 1 sel ext(T, ψ, tb, t, !TBV )@ .P || s2 ::ρ 2 (I, R) → s1 ::ρ 1 P σ || s2 ::ρ 2 (I, R) (SEL INT) I = I[id → ⊥][sk → I.sk ↓ T t ] σ = [(I , R )/TBV ] R = {E tσ ρ 1 | ∃t : t ∈ R ∧ match(E T ρ 1 , t ) ∧ σ = [t /E T ρ 1 ] ∧ ψσ} s1 ::ρ 1 sel int(T, ψ, (I, R), t, !TBV ).P → s1 ::ρ 1 P σ (UPD) ρ1(l) = s2 I.id = tb R 1 = {t | t ∈ R ∧ ¬(∃σ : match(E T ρ 1 , t ) ∧ σ = [t /E T ρ 1 ] ∧ ψσ ∧ E tσ ρ 1 |= I.sk)} R 2 = {E tσ ρ 1 | ∃t : t ∈ R ∧ match(E T ρ 1 , t ) ∧ σ = [t /E T ρ 1 ] ∧ ψσ ∧ E tσ ρ 1 |= I.sk}
s1 ::ρ 1 update(T, ψ, t, tb)@ .P || s2 ::ρ 2 (I, R) → s1 ::

ρ 1 P ||s2 ::ρ 2 (I, R 1 R 2 ) (AGGR) t = f ({t | ∃σ : t ∈ R ∧ match(E T ρ 1 , t) ∧ σ = [t/E T ρ 1 ] ∧ ψσ ) ρ(l) = s2 I.id = tb match(E T ρ 1 , t )
s1 ::ρ 1 aggr(T, ψ, tb, f, T )@ .P || s2 ::ρ 2 (I, R) → s1 :: We proceed with a detailed explanation of the semantic rules in Table 3 that account for the execution of Klaim-DB actions and the ones in Table 4 that mainly describe the execution of the control flow constructs. In the explanation, we will avoid reiterating that each table resides in a database located at some , but directly state "table ... located at ".

ρ 1 P [t /E T ρ 1 ] || s2 ::ρ 2 (I, R) (CREATE) ρ1 (

Semantics for Actions

Insertion and Deletion The rule (INS) of Table 3 says that to insert a row t into a table tb at s 2 , the logical locality needs to be evaluated to s 2 under the local environment ρ 1 , the table identifier tb must agree with that of a destination table (I, R) already existing at s 2 , and the tuple t needs to satisfy the requirements imposed by the schema I.sk. If these conditions are met, then the evaluated tuple is added into the data set R.

Example 2 (Adding new Shoes). The local action at s 1 that inserts an entry for KLD high boots of a new color, white, sized "37", produced in 2015, with 6 in stock, is insert(T 0 , KLD)@self , where T 0 = ("001", "HighBoot", "2015", "white", "37", 6, 0). By (INS), we have

s1 ::ρ 1 insert(T0, KLD)@self .nil || s1 ::ρ 1 (I1, R1) → s1 ::ρ 1 nil || s1 ::ρ 1 (I1, R 1 ), where R 1 = R 1 {T 0 }, since ρ 1 (self ) = s 1 , KLD = I 1 .id and T 0 |= I 1 .sk .
Deletion operations are performed according to the rule (DEL). A deletion can be performed if the logical locality refers to the physical locality s 2 under the local environment ρ 1 , and the specified table identifier tb agrees with that of the table (I, R) targeted. The rows that do not match the pattern T or do not satisfy the condition ψ will constitute the resulting data set of the target.

Example 3 (Deleting Existing Shoes). The local action at 1 , deleting all entries for white KLD high boots sized "37" produced in 2015, from the resulting table of Example 2, is delete(T 0 , true, KLD)@self , where T 0 = ("001", "HighBoot", "2015", "W hite", "37", !x, !y).

We have the transition:

s1 ::ρ 1 delete(T0, true, KLD)@self .nil || s1 ::ρ 1 (I1, R 1 ) → s1 ::ρ 1 nil || s1 ::ρ 1 (I1, R1).
This reflects that the original table is recovered after the deletion.

Selection, Update and Aggregation The rule (SEL EXT) describes the way the "external" selection operations are performed. It needs to be guaranteed that the logical locality is evaluated to s 2 (under the local environment ρ 1 ), and the identifier tb is identical to that of the table (I, R) existing at s 2 . If the conditions are met, all the rows that match the pattern T and satisfy the predicate ψ, will be put into the result data set R . The schema of the resulting table is also updated according to the pattern T and the tuple t. The resulting table (I , R ) is substituted into each occurrence of the table variable TBV used in the continuation P .

Example 4 (Selection of Shoes in a Certain Color). The Klaim-DB action performed from the head office, selecting the color, size, and sales of the types of high boots that are not red, at the local database at s 1 , is sel ext(("001", "HighBoot", !x, !y, !z, !w, !p), y = "red", KLD, (y, z, p), !TBV )@ 1.

According to the rule (SEL EXT), we have the transition:

s0 ::ρ 0 sel ext(("001", "HighBoot", !x, !y, !z, !w, !p), y = "red", KLD, (y, z, p), !TBV )@ 1.nil || s1 ::ρ 1 (I1, R1) → s0 ::ρ 0 nil || s1 ::ρ 1 (I1, R1).

The conditions ρ 0 ( 1 ) = s 1 and KLD = I.id in the premises of the rule are satisfied. The I in (SEL EXT) is such that I .id = ⊥, I .sk =< "Color" : String, "Size" : Int, "Sales" : Int >, and I agrees with I on the other attributes. The table variable TBV is replaced by (I , R ), for some R = {("black", "38", "2"), ("black", "37", "2")}.

"Internal" selections are performed according to the rule (SEL INT). No constraints concerning localities/identifiers are needed. All the rows that match the pattern T and satisfy the predicate ψ are selected into the resulting data set R . The schema of the result table is produced by using the projection operation introduced in the beginning of this section. All occurrences of the table variable TBV in the continuation P will be replaced by the resulting table (I , R ) before continuing with further transitions. We use ⊥ for I .id to indicate that (I , R ) is a "temporary" table that can be thought of as existing in the query engine, rather than in the database.

The updates of data stored in tables are executed according to the rule (UPD). Again, it is checked that the specified logical locality of the table (I, R) to be updated corresponds to the site where the table actually resides, and that the specified table identifier matches the actual one. An update goes through all the elements t of R. This t is updated only if it matches the template T , resulting in the substitution σ that makes the predicate ψ satisfied, and tσ is evaluated to a tuple that satisfies the schema I.sk. The evaluation result of tσ will then replace t in the original table, which is captured by R 2 . If at least one of the above mentioned conditions is not met, then t will be left intact (described by R 1 ).

Example 5 (Update of Shoes Information). Suppose two more red KLD high boots sized 37 are sold. The following local action informs the database at s 1 of this. update( ("001", "HighBoot", "2015", "red", "37", !x, !y), true, ("001", "HighBoot", "2015", "red", "37", x -2, y + 2), KLD)@self .

We have the transition:

s1 ::ρ 1 update(("001", "HighBoot", "2015", "red", "37", !x, !y), true, ("001", "HighBoot", "2015", "red", "37", x -2, y + 2), KLD)@self .nil || s1 ::

ρ 1 (I, R) → s1 ::ρ 1 nil || s1 ::ρ 1 (I, R 1 R 2 ).
The multiset R 1 consists of all the entries that are intact -shoes that are not red high boots sized 37, while R 2 contains all the updated items.

The rule (AGGR) describes the way aggregation operations are performed. The matching of localities and table identifiers is still required. The aggregation function f is applied to the multiset of all tuples matching the template T , as well as satisfying the predicate ψ. The aggregation result t obtained by the application of f is bound to the specified template T only if the evaluation of T matches t . In that case the substitution reflecting the binding is applied to the continuation P and the aggregation is completed.

Example 6 (Aggregation of Sales Figures).

The local aggregation returning the total sales of the shoes with ID "001", can be modelled by the Klaim-DB action aggr(T 0 , true, KLD, sum 7 , (!res))@self , where T 0 = ("001", !x, !y, !z, !w, !q, !o), and sum 7 = λR.(sum({v 7 |(v 1 , ..., v 7 ) ∈ R})), i.e., sum 7 is a function from multisets R to unary tuples containing the summation results of the 7-th components of the tuples in R.

It is not difficult to derive:

s1 ::ρ 1 aggr(T0, true, KLD, sum7, (!res))@ 1.nil || s1 ::ρ 1 (I, R) → s1 ::ρ 1 nil || s1 ::ρ 1 (I, R).

The variable !res is then bound to the integer value 12 (2 + 5 + 1 + 2 + 2).

Example 7 (Selection using Aggregation Results

). Consider the query from s 0 that selects all the colors, sizes and sales of the types of high boots produced in the year 2015, whose ID is "001" and whose sales figures are above average. This query can be modelled as a sequence of actions at s 0 , as follows.

s0 ::ρ 0 aggr(T0, true, KLD, avg 7 , !res)@ 1.sel ext(T0, w ≥ res, KLD, (x, y, w), !TBV )@ 1.nil || s1 ::ρ 1 (I1, R1), where T 0 = ("001 , "HighBoot", "2015", !x, !y, !z, !w) and avg 7 = λR.( sum({v7|(v1,...,v7)∈R})

|R|

) is a function from multisets R to unary tuples containing the average value of the 7-th components of the tuples in R.

Adding and Dropping Tables To create a new table, the rule (CREATE) ensures that the physical site s 2 corresponding to the logical locality mentioned in the action create(I)@ does exist. Then a table with the interface I and an empty data set is created at the specified site. It remains an issue to ascertain that there are no existing tables having the same identifier I.id at the target site. This is achieved with the help of the rule (PAR) in Table 4.

To drop an existing table, the rule (DROP) checks that a table with the specified identifier tb does exist at the specified site ρ 1 (l) = s 2 . Then the table is dropped by replacement of it with nil.

Semantics for Processes and Networks

Directing our attention to Table 4, the rules (FOR tt s ) and (FOR ff s ) specify when and how a loop is to be executed one more time, and has finished, respectively.

(FOR tt s ) t0 ∈ R match(E T ρ, t0) s ::ρ foreachs T in R : P → s ::ρ P [t0/E T ρ]; foreachs T in R \ {t0} : P (FOR tt p ) t0 ∈ R match(E T ρ, t0) s ::ρ foreachp T in R : P → s ::ρ P [t0/E T ρ] | foreachp T in R \ {t0} : P (FOR ff s ) ¬(∃t0 ∈ R : match(E T ρ, t0)) s ::ρ foreachs T in R : P → s ::ρ nil (FOR ff p ) ¬(∃t0 ∈ R : match(E T ρ, t0)) s ::ρ foreachp T in R : P → s ::ρ nil (SEQ tt ) s ::ρ P1 || N → s ::ρ P 1 || N s ::ρ P1; P2 || N → s ::ρ P 1 ; P2 || N (SEQ ff ) s ::ρ P1 || N → s ::ρ nil || N s ::ρ P1; P2 || N → s ::ρ P2 || N (CALL) s ::ρ A( t) → s ::ρ P [ṽ/x] if A(x) P ∧ E t ρ = ṽ (PAR) N1 → N 1 N1||N2 → N 1 ||N2 if Lid (N 1 ) ∩ Lid (N2) = ∅ (RES) N → N (νs)N → (νs)N (EQUIV) N1 ≡ N2 N2 → N3 N3 ≡ N4 N1 → N4
The rule (FOR tt s ) says that if there are still more tuples (e.g., t 0 ) in the multiset R matching the pattern specified by the template T , then we first execute one round of the loop with the instantiation of variables in T by corresponding values in t 0 , and then continue with the remaining rounds of the loop by removing t 0 from R. The rule (FOR ff s ) says that if there are no more tuples in R matching T , then the loop is completed. The rules (FOR tt p ) and (FOR ff p ) can be understood similarly.

The rules (SEQ tt ) and (SEQ ff ) describe the transitions that can be performed by sequential compositions P 1 ; P 2 . In particular, the rule (SEQ tt ) accounts for the case where P 1 cannot finish after one more step; whereas the rule (SEQ ff ) accounts for the opposite case. We require that no variable bound in P 1 is used in P 2 , thereby avoiding the need for recording the substitutions resulting from the next step of P 1 , that need to be applied on P 2 .

The rule (PAR) says that if a net N 1 can make a transition, then its parallel composition with another net N 2 can also make a transition, as long as no clashes of table identifiers local to each site in N 1 ||N 2 are introduced. It makes use of the function Lid (...) that gives the multiset of pairs of physical localities and table identifiers in networks and components. This function is overloaded on networks and components, and is defined inductively as follows.

Lid (nil) = ∅ Lid (s, P ) = ∅ Lid (N1||N2) = Lid (N1) Lid (N2) Lid (s, (I, R)) = {(s, I.id)} Lid ((νs)N ) = Lid (N ) Lid (s, C1|C2) = Lid (s, C1) Lid (s, C2) Lid (s ::ρ C) = Lid (s, C)
The rules (CALL), (RES) and (EQUIV) are self-explanatory. It is a property of our semantics that no local repetition of table identifiers can be caused by a transition (including the creation of a new table). Define no rep(A) = (A = set(A)), which expresses that there is no repetition in the multiset A, thus A coincides with its underlying set. This property is formalized in Lemma 1 whose proof is straightforward by induction on the derivation of the transition.

Lemma 1. For all nets N and N , if no rep(Lid (N )) and N → N , then it holds that no rep(Lid (N )).

Thus, imposing "non-existence of local repetition of table identifiers" as an integrity condition for the initial network will guarantee the satisfaction of this condition in all possible derivatives of the network. Suppose there is no other table at s 1 with the identifier KLD. By (PAR), and the structural congruence, we have where t 1 = ("001", "HighBoot", "2015", "white", "37", "6", "0").

Case Study

Continuing with our running example, we illustrate the modelling of data aggregation over multiple databases local to different branches of the department store chain in Klaim-DB. In more detail, a manager of the head office wants statistics on the total sales of KLD high boots from the year 2015, in each branch operating in Copenhagen.

We will think of a procedure stat at the site s 0 of the head office, carrying out the aggregation needed. Thus the net for the database systems of the department store chain, as considered in Section 2, specializes to the following, where C 0 is the remaining tables and processes at s 0 apart from Stores and stat. A detailed specification of the procedure stat is then given in Figure 4. First of all, a result table with the interface I res is created, where I res .id = result and I.sk =< "Brand" : String, "City" : String, "Shop name" : String, stat create(Ires)@self .

sel ext((!x, !y, !z, !w, !p), KLD ∈ w ∧ x = "CP H", Stores, (z, p), !TBV )@self . foreachp (!q, !u) in TBV .R : aggr(("001", "HighBoot", "2015", !x, !y, !z, !w), true, KLD, sum7, (!res))@u. insert((q, "HighBoot", res), result)@self ; ... drop(result)@self Fig. 4. The Procedure for Distributed Data Aggregation "Sales" : Int >. Then all the logical localities of the local databases used by the branches in Copenhagen that actually sell KLD shoes are selected, together with the shop names of such branches. This result set is then processed by a parallel loop. The number of KLD high boots from 2015 that are sold is counted at each of these localities (branches), and is inserted into the resulting table together with the corresponding shop name and the information "High Boot" describing the shoe type concerned. The resulting table, displayed in Figure 5, can still be queried/manipulated before being dropped. In the end of this case study, we would like to remark that the use of the parallel loop in carrying out all the individual remote aggregations has made the overall query more efficient. Hence some performance issues can be captured by modelling in Klaim-DB.

Shop name

Extension and Discussion

Joining Tables Our modelling language can be extended to support querying from the "join" of tables. For this extension, we make use of under-specified join expressions je, that can have table identifiers or table variables as parameters. The syntax of the new selection actions is shown below.

a ::= ... | sel ext(T, ψ, je(tb1, ..., tbn), t, !TBV )@ 1, ..., n | sel int(T, ψ, je(TBV 1, ..., TBV n), t, !TBV ) | ...

For an external selection, we allow to use a list 1 , ..., n of localities (abbreviated as ), to join tables located in multiple databases. In more detail, for i ∈ {1, ..., n}, tb i is supposed to be the identifier of a table located at i . We use je R and je I to represent the interpretation of je in terms of how the datasets and the schemas of the joined tables should be combined. As an example, a plain list of table identifiers or tables (substitution of table variables) corresponds to taking the concatenation of all the schemas and selecting from the cartesian product of all the data sets.

n = |tb| = | | ∧ ∀j ∈ {1, ..., n} : ρ0(lj) = sj ∧ tbj = Ij.id I = [id → ⊥][sk → ( je I (I)) ↓ T t ] σ = [(I , R )/TBV ] R = {E tσ ρ 1 | ∃t : t ∈ je R (R) ∧ match(E T ρ 1 , t ) ∧ σ = (t /E T ρ 1 ) ∧ ψσ}
(SELJ INT) I = [id → ⊥][sk → ( je I (I)) ↓ T t ] σ = [(I , R )/TBV ] R = {E tσ ρ 1 | ∃t : t ∈ je R (R) ∧ match(E T ρ 1 , t ) ∧ σ = (t /E T ρ 1 ) ∧ ψσ} s1 ::ρ 1 sel int(T, ψ, je((I, R)), t, !TBV ).P → s1 ::ρ 1 P σ
The pattern matching against the template T and the satisfaction of the predicate ψ will be examined on tuples from je R (R 1 , ..., R n ), and the predicate ψ can now impose constraints on the fields from different tables.

The adaption needed for the semantics is fairly straightforward. The semantic rules (SELJ EXT) and (SELJ INT) that describe selection operations from joined tables are presented in Table 5. In the rule (SELJ EXT), although it is stipulated that the j-th table identifier specified in the list tb must be identical to the j-th table (I j , R j ) listed as parallel components, no undesired stuck configurations are caused because of the structrual congruence.

This extension paves the way for the general ability to operate on multiple databases by a single action, which is in line with the design philosophy of multidatabase systems (e.g., [START_REF] Kuhn | Vip-mdbs: A logic multidatabase system[END_REF]).

Discussion

We could have required whole tables to be received into local variables by using the standard Klaim actions out(t)@ and in(t)@ , and made the selection and aggregation operations work only on table variables. In this way we could have gotten rid of "external selection". However, "external selection" on remote localities can potentially reduce the communication cost considerably, since only one tuple (for aggregation) or a part of a table (for selection) need to be returned. For selection the reduction is particularly meaningful when the resulting data set is small.

Concerning the result of selection operations, an alternative that we have not adopted is the direct placement of the result in a separate table. This table is either created automatically by the selection operation itself, with an identifier specified in the selection action, or a designated "result table" at each site. However, a problem is that the removal of the automatically created tables will need to be taken care of by the system designer making the specification, using drop(I)@ actions. And similar problems arise with the maintenance of the designated "result table" local to each site (e.g., the alteration of its schema, the cleaning of old results, etc.). To abstain from these low-level considerations, table variables are introduced and binding is used for the selection results.

The interoperability between database systems and ordinary applications can also be realized by bringing back the primitive Klaim actions (out(t)@ , in(T )@ , and eval(P )@ ) and allowing the co-existence of tables and plain tuples at different localities. Via the re-introduction of the eval action, we would also be able to send out mobile processes to perform complex data processing on-site.

Conclusion

We have proposed Klaim-DB -a modelling language that borrows insights from both the coordination language Klaim and database query languages like SQL, to support the high-level modelling and reasoning of distributed database applications. The semantics is illustrated with a running example on the query and management of the databases used by a department store chain. Data aggregation across the geographically scattered databases at the individual stores, as performed by a coordinator, is then modelled in the language. In the model, the local aggregations at the store-owned databases are performed in parallel with each other, benefiting the performance.

Our use of templates in the query actions of Klaim-DB is in line with the spirit of the QBE language (Query by Example [START_REF] Moshé | Query by example[END_REF]). The choice of using multiset operations for the semantics of these actions, on the other hand, has the flavor of the Domain Relational Calculus underlying QBE. The work presented in [START_REF] Madhan | On the use of linda as a framework for distributed database systems[END_REF] discusses typical coordination issues in distributed databases in Linda, from an informal, architectural viewpoint. This work is marginally related to ours.

The specification and enforcement of security policies is an important concern in distributed database systems. A simple example is: when inserting data into a remote table, it is important to know whether the current site trusts the remote site with respect to confidentiality, and whether the remote site trusts the current site with respect to integrity. Recently, [START_REF] Dwork | Differential privacy[END_REF] and [START_REF] Lourenço | Information flow analysis for valued-indexed data security compartments[END_REF] address privacy and information flow issues in database applications. Another work, [START_REF] Terkel | Locality-based security policies[END_REF], provides an information flow analysis for locality-based security policies in Klaim. By elaborating on the language design, we provide in this paper a solid ground for any future work aiming to support security policies and mechanisms.

Another interesting line of future work is the specification of transactions (e.g., [START_REF] Kühn | A parallel logic language for transaction specification in multidatabase systems[END_REF]). This is needed for the modelling of finer-grained coordination between database accesses.

Fig. 1 .

 1 Fig. 1. Running Example

  l) = s2 s1 ::ρ 1 create(I)@ .P ||s2 ::ρ 2 nil → s1 ::ρ 1 P ||s2 ::ρ 2 (I, ∅) (DROP) ρ1(l) = s2 I.id = tb s1 ::ρ 1 drop(tb)@ .P ||s2 ::ρ 2 (I, R) → s1 ::ρ 1 P ||s2 ::ρ 2 nil

Example 8 (

 8 Transition of Networks). Continuing with Example 2, we illustrate how our semantics help establish global transitions.

  s0 ::ρ 0 (I0, R0)|C 0 || s1 ::ρ 1 (I1, R1)|insert(t1, KLD)@self .nil|C 1 || ... || sn ::ρ n (In, Rn)|C n → s0 ::ρ 0 (I0, R0)|C 0 || s1 ::ρ 1 (I1, R 1 )|C 1 || ... || sn ::ρ n (In, Rn)|C n ,

  s0 ::ρ 0 (stat|(I, R)|C 0 ) || s1 ::ρ 1 ((I1, R1)|C 1 ) || ... || sn ::ρ n ((In, Rn)|C n )

Fig. 5 .

 5 Fig. 5. The Table result

  je I (I1, ..., In) = j∈{1,...,n} Ij.sk je R (R1, ..., Rn) = R1 × ... × Rn

  s0 ::ρ 0 sel ext(T, ψ, je(tb), t, !TBV )@ .P || s1 ::ρ 1 (I1, R1) || ... || sn ::ρ n (In, Rn) → s0 ::ρ 0 P σ || s1 ::ρ 1 (I1, R1) || ... || sn ::ρ n (In, Rn)

Table 1 .

 1 The Syntax of Klaim-DB

	N

  as the local environment for each site s j -this restricts database accesses across different local databases and corresponds to a centralized architecture.

	City	Address	Shop Name	Brand	Logical Locality
	CPH	ABC DEF 2, 1050	Shop1	{KLD, SH, ...}	1
	CPH	DEF HIJ 13, 2800	Shop2	{KLD, SH, ...}	2
	CPH HIJ KLM 26, 1750	Shop3	{KLD, SH, ...}	3
	AAL KLM NOP 3, 3570	Shop4	{LAM, IM K, ...}	4
	AAL NOP QUW 18, 4500	Shop5	{LAM, IM K, ...}	5
	...	...	...	...	...

Table 3 .

 3 The Semantics for Actions

	(INS)

Table 4 .

 4 The Semantics for Processes and Nets (Continued)

Table 5 .

 5 The Semantic Rules for Selection from Joins

	(SELJ EXT)

Acknowledgment. We would like to thank the anonymous reviewers for their valuable feedback. Xi Wu is partly supported by the IDEA4CPS project, ECNU Project of Funding Overseas Short-term Studies, Domestic Academic Visit and NSFC Project (No.61361136002 and No.61321064). Ximeng Li is partly funded by SESAMO, a European Artemis project.