domain: info math

text:
S Graillat
email: stef.graillat@upmc.fr

F Jézéquel

R Picot

Numerical validation of compensated algorithms with stochastic arithmetic

Keywords: CADNA, compensated algorithms, Discrete Stochastic Arithmetic, error-free transformations, floating-point arithmetic, numerical validation, rounding errors

Compensated algorithms consist in computing the rounding error of individual operations and then adding them later on to the computed result. This makes it possible to increase the accuracy of the computed result efficiently. Computing the rounding error of an individual operation is possible through the use of a so-called error-free transformation. In this article, we show that it is possible to validate the result of compensated algorithms using stochastic arithmetic. We study compensated algorithms for summation, dot product and polynomial evaluation. We prove that the use of the random rounding mode inherent to stochastic arithmetic does not change the accuracy of compensated methods. This is due to the fact that error-free transformations are no more exact but still sufficiently accurate to improve the numerical quality of results.

Introduction

Computing power rapidly increases and Exascale computing (10 18 floating-point operations per second) should be reached in a few years. Such a computing power also means a lot of rounding errors. Indeed, nearly all floating-point operations imply a small rounding which can accumulate along the computation and finally an incorrect result may be produced. As a consequence, it is fondamental to be able to give some information about the numerical quality of the computed results. By numerical quality, we mean here the number of significant digits of the computed result that are not affected by rounding errors.

Definitions and notations

In this paper, we assume to work with a binary floating-point arithmetic adhering to IEEE 754 floating-point standard [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF] and we suppose that no overflow occurs. The set of floating-point numbers is denoted by F, the relative rounding error by u. For IEEE 754 double precision, we have u = 2 -53 and for single precision u = 2 -24 .

We denote by fl * (•) the result of a floating-point computation, where all operations inside parentheses are done in floating-point working precision with a directed rounding (that is to say toward -∞ or +∞). Floating-point operations in IEEE 754 satisfy [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF] ∃ We use standard notations for error estimations. The quantities γ n are defined as usual [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF] by

ε 1 ∈ R, ε 2 ∈ R such that fl * (a • b) = (a • b)(1 + ε 1 ) = (a • b)/(1 + ε 2 ) for • = {+, -} and |ε ν | ≤ 2u. ( 2
γ n (u) := nu 1 -nu for n ∈ N,
where implicitly assumed that nu < 1.

To keep track of the (1 + ε) factors in our error analysis, we use the relative error counters introduced by Stewart [START_REF] Stewart | Introduction to matrix computations[END_REF]. For a positive integer n, n denotes the following product

n (u) = n i=1
(1 + ε i ) ρ i with ρ i = ±1 and |ε i | ≤ u (i = 1, . . . , n).

The relative error counters verify j (u) k (u) = j (u)/ k (u) = j + k (u). When n denotes any error counter, then there exists a quantity θ n such that n (u) = 1 + θ n (u) and |θ n (u)| ≤ γ n (u).

Remark 1. We give the following relations about γ n , that will be frequently used in the sequel of the paper. For any positive integer n, nu ≤ γ n (u), γ n (u) ≤ γ n+1 (u), (1 + u)γ n (u) ≤ γ n+1 (u), 2nu(1 + γ 2n-2 (u)) ≤ γ 2n (u).

Principles of Discrete Stochastic Arithmetic (DSA)

This section briefly recalls the principles of Discrete Stochastic Arithmetic. Further details are given for instance in [START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF]. Based on the CESTAC method [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF], Discrete Stochastic Arithmetic enables one to estimate in a computed result which digits are affected by rounding errors. It requires a random rounding mode that consists in rounding any result upwards or downwards with the same probability.

To use the CESTAC method in a code that computes a result R, one executes N times this code with the random rounding mode. Therefore N different results R i are obtained. The value of the computed result R is chosen to be the mean value of {R i } and, if no overflow occurs, the number of exact significant digits in R can be estimated as

C R = log 10 √ N R στ β
where σ is the standard deviation of {R i } and τ β is the value of Student's distribution for N -1 degrees of freedom and a confidence level 1 -β. In practice β = 0.05 and N = 3. Indeed, it has been shown [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF][START_REF] Chesneaux | Sur la robustesse de la méthode CESTAC[END_REF] that N = 3 is in some reasonable sense the optimal value. The estimation with N = 3 is more reliable than with N = 2 and increasing the size of the sample does not improve the quality of the estimation.

The CESTAC method relies on a first order model of rounding errors which becomes invalid if both operands in a multiplication or the divisor in a division are not significant [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF]. Therefore the CESTAC method requires to control all multiplications and divisions during the execution. This so-called self-validation has led to the concept of computational zero [START_REF] Vignes | Zéro mathématique et zéro informatique[END_REF], defined below, and also to the synchronous implementation of the method: each arithmetic operation is performed N times before the next one is executed. Definition 3.1. A result R = {R i } computed using the CESTAC method is a computational zero, denoted by @.0, if ∀i, R i = 0 or C R ≤ 0. A computational zero is either zero or a result that has no more correct digits because of rounding errors. From the concept of computational zero, discrete stochastic relations [START_REF] Chesneaux | Les fondements de l'arithmétique stochastique[END_REF] have been defined as follows.

Definition 3.2. Let X = {X i } and Y = {Y i } be two results computed with the CESTAC method, 1. X = Y if and only if X -Y = @.0, 2. X > Y if and only if X > Y and X -Y = @.0, 3. X ≥ Y if and only if X ≥ Y or X -Y = @.0. Discrete Stochastic Arithmetic (DSA) is the combination of the CESTAC method, the concept of computational zero, and the discrete stochastic relations [START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF].

The CADNA 1 (Control of Accuracy and Debugging for Numerical Applications) library [START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF] is an implementation of DSA devoted to programs written in C/C++ and Fortran. CADNA provides, with the probability 1 -β = 95%, the number of exact significant digits of any computed result. The CADNA library allows to use new numerical types: the stochastic types. Each stochastic variable contains N = 3 values of the corresponding floating-point type, one for each sample R i . When a stochastic variable is printed, only its exact significant digits appear. Arithmetic operators, comparison operators, all the mathematical functions have been overloaded to return a stochastic type when called with stochastic arguments. Therefore the use of CADNA in a program requires only a few modifications: essentially changes in the declarations of variables and in input/output statements. During the execution, CADNA can detect numerical instabilities, which are usually due to the presence of numerical noise. When numerical instabilities are detected, dedicated CADNA counters are incremented. At the end of the run, the value of these counters together with appropriate warning messages are printed on standard output.

In the next sections, we aim at analysing the effects of the random rounding mode required by DSA on compensated algorithms a priori intended to be used with rounding to nearest. Therefore we will present the impact of a directed rounding mode on the accuracy of results provided by compensated algorithms.

Accurate summation

This section presents the accuracy obtained with the classic summation algorithm and with various compensated summation algorithms, using either rounding to nearest or directed rounding.

Classic summation

The classic algorithm for computing summation is the recursive Algorithm 1.

function res = Sum(p) 1: s 1 ← p 1 2: for i = 2 to n do 3: s i ← s i-1 + p i 4: end for 5: res ← s n Algorithm 1: Summation of n floating-point numbers p = {p i }
The error generated by Algorithm 1 is recalled in Proposition 4.1.

|res -s| |s| ≤ γ n-1 (u) cond p i .
With directed rounding, if nu < 1 2 , the result res of Algorithm 1 satisfies

|res -s| |s| ≤ γ n-1 (2u) cond p i .
Because γ n-1 (u) ≈ (n-1)u, the bound for the relative error is essentially nu times the condition number. This accuracy is sometimes not sufficient in practice. Indeed, when the condition number is large (greater than 1/u) then the recursive algorithm does not even return one correct digit. Algorithms to evaluate more accurately the sum of floating-point numbers are presented in the sequel of this section.

Compensated summation with rounding to nearest

Error-free transformations exist for the sum of two floating-point numbers with rounding to nearest: TwoSum [START_REF] Knuth | The Art of Computer Programming[END_REF] which requires 6 floating-point operations and FastTwoSum [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], given as Algorithm 2, which requires a test and 3 floating-point operations. These algorithms compute both the floatingpoint sum c of two numbers a and b and the associated rounding error d such that c+d = a+b. Another algorithm, proposed by Priest in [21, p.14-15] and given later as Algorithm 4 (PriestTwoSum), although more costly, computes with any rounding mode an error-free transformation for the sum of two floating-point numbers.

A compensated algorithm to evaluate accurately the sum of n floating-point numbers is presented as Algorithm 3 (FastCompSum) [START_REF] Neumaier | Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen[END_REF]. This sum is corrected thanks to an error-free transformation used for each individual summation. Although FastTwoSum is called in Algorithm 3, with rounding to nearest the same result can be obtained using another error-free transformation (TwoSum or PriestTwoSum).

The error on the result res of Algorithm 3 obtained with rounding to nearest is analysed in [START_REF] Ogita | Accurate sum and dot product[END_REF]. A bound for the absolute error is recalled in Proposition 4.

1: π 1 ← p 1 2: σ 1 ← 0 3: for i = 2 to n do 4: [π i , q i ] ← FastTwoSum(π i-1 , p i ) 5:
σ i ← σ i-1 + q i 6: end for [START_REF] Ogita | Accurate sum and dot product[END_REF]). With rounding to nearest, if nu < 1, then, also in the presence of underflow, the result res of Algorithm 3 (FastCompSum) satisfies

|res -s| |s| ≤ u + γ 2 n-1 (u) cond p i .
From Corollary 4.4, because γ n-1 (u) ≈ (n -1)u, the bound for the relative error on the result is essentially (nu) 2 times the condition number plus the rounding u due to the working precision. The second term on the right hand side reflects that the computation is carried out as in twice the working precision (u 2 ). The first term represents the rounding back in the working precision.

Compensated summation with directed rounding

We recall here the impact of a directed rounding mode on Algorithm 3 (FastCompSum). With directed rounding, Algorithm 2 (FastTwoSum) is not an error-free transformation. A bound on the difference between the floating-point number d computed by Algorithm 2 and the error e due to the floating-point addition is recalled in Proposition 4.5. function res = PriestCompSum(p)

end if 4: c ← a + b 5: e ← c -a 6: g ← c -e 7: h ← g -a 8: f ← b -h 9: d ← f -e 10: if d + e = f
1: π 1 ← p 1 2: σ 1 ← 0 3: for i = 2 to n do 4: [π i , q i ] ← PriestTwoSum(π i-1 , p i ) 5:
σ i ← σ i-1 + q i 6: end for 7: res ← π n + σ n Algorithm 5: Compensated summation of n floating-point numbers p = {p i } using PriestTwoSum Lemma 4.9 is given for later use in Sect. 7. Lemma 4.9. Let us suppose Algorithm 5 (PriestCompSum) is applied, with directed rounding, to floating-point numbers

p i ∈ F, 1 ≤ i ≤ n. Let s := p i and S := |p i |. If nu < 1 2 , then n i=2 |q i | + |π n | ≤ |s| + 2γ n-1 (2u)S. (4.6)
Proof.

As

n i=2 |q i | + |π n | = n i=2 |q i | + |s - n i=2 q i |, we have n i=2 |q i | + |π n | ≤ |s| + 2 n i=2 |q i |. (4.7)
From Lemma 4.6, we deduce that

n i=2 |q i | ≤ γ n-1 (2u) n i=1 |p i | = γ n-1 (2u)S. (4.8)
Finally Equation 4.6 is obtained from Equations 4.7 and 4.8.

A bound for the absolute error on the result of Algorithm 5 (PriestCompSum) obtained with directed rounding is given in Proposition 4.10. Proposition 4.10. Let us suppose Algorithm PriestCompSum is applied, with directed rounding, to floating-point numbers

p i ∈ F, 1 ≤ i ≤ n. Let s := p i and S := |p i |. If nu < 1 2 , then, also in the presence of underflow, |res -s| ≤ 2u|s| + γ 2 n-1 (2u)S. (4.9)
The proof is similar to the one given in [START_REF] Ogita | Accurate sum and dot product[END_REF] for compensated summation with rounding to nearest.

Proof. Because σ n = fl * ( n i=2 q i ), we have

|σ n - n i=2 q i | ≤ γ n-2 (2u) n i=2 |q i |.
Therefore, from Equation 4.8, we deduce that

|σ n - n i=2 q i | ≤ γ n-2 (2u)γ n-1 (2u)S.
As Algorithm 5 is executed with directed rounding, it yields

res = fl * (π n + σ n ) = (1 + ε)(π n + σ n ) with |ε| ≤ 2u, |res -s| = | fl * (π n + σ n ) -s|, |res -s| = |(1 + ε)(π n + σ n -s) + εs|, |res -s| = |(1 + ε)(π n + n i=2 q i -s) + (1 + ε)(σ n - n i=2 q i ) + εs|. Since s = n i=1 p i = π n + n i=2 q i , (4.10)
then

|res -s| ≤ (1 + 2u)|σ n - n i=2 q i | + 2u|s|, |res -s| ≤ (1 + 2u)γ n-2 (2u)γ n-1 (2u)S + 2u|s|. ( 4.11)
We have

(1 + 2u)γ n-1 (2u) < γ n (2u). (4.12)
Indeed, it is clear that

γ n (2u) -(1 + 2u)γ n-1 (2u) = 2u 1 + 2nu (1 -2(n -1)u) (1 -2nu) ,
and then

γ n (2u) -(1 + 2u)γ n-1 (2u) > 0.
Finally Equation 4.9 can be deduced from Equations 4.11 and 4.12.

From Proposition 4.10, a bound for the relative error on the result of Algorithm 5 (PriestCompSum) obtained with directed rounding is deduced in Corollary 4.11.

Corollary 4.11. With directed rounding, if nu < 1 2 , then, also in the presence of underflow, the result res of Algorithm 5 (PriestCompSum) satisfies

|res -s| |s| ≤ 2u + γ 2 n-1 (2u) cond p i .
Like with Algorithm 3 (FastCompSum), we deduce from Corollary 4.11 that the relative error bound on the result of Algorithm 5 (PriestCompSum) computed with directed rounding is essentially (nu) 2 times the condition number plus the rounding 2u due to the working precision.

Accurate dot product

In this section, we present the accuracy obtained with the classic dot product algorithm. We also present an algorithm that enables one to compute a dot product as in twice the working precision with rounding to nearest [START_REF] Ogita | Accurate sum and dot product[END_REF]. We recall the error on its result computed with rounding to nearest. Then we analyse the impact of a directed rounding mode on this algorithm. In this section, we assume that no underflow occurs.

Classic dot product

The classic algorithm for computing a dot product is Algorithm 6. function res = Dot(x, y)

1: s 1 ← x 1 y 1 2: for i = 2 : n do 3: s i ← x i • y i + s i-1 4: end for 5: res ← s n Algorithm 6: Classic dot product of x = {x i } and y = {y i }, 1 ≤ i ≤ n
The following proposition sums up the properties of this algorithm.

Proposition 5.1. Let floating point numbers x i , y i ∈ F, 1 ≤ i ≤ n, be given and denote by res ∈ F the result computed by Algorithm 6 (Dot). With rounding to nearest, if nu < 1, we have

|res -x T y| ≤ γ n (u)|x T ||y|, (5.13)
and with directed rounding, if nu < 1 2 , we have

|res -x T y| ≤ γ n (2u)|x T ||y|. (5.14)
Proof. The proof can be found in Higham [11, p.63].

We can rewrite the previous inequalities in terms of the condition number of the dot product defined by cond(x T y) = 2 |x| T |y| |x T y| .

Corollary 5.2. With rounding to nearest, if nu < 1, the result res of Algorithm 6 satisfies

|res -x T y| |x T y| ≤ 1 2 γ n (u) cond(x T y).
With directed rounding, if nu < 1 2 , the result res of Algorithm 6 satisfies

|res -x T y| |x T y| ≤ 1 2 γ n (2u) cond(x T y).

Compensated dot product with rounding to nearest

A compensated dot product algorithm is presented in [START_REF] Ogita | Accurate sum and dot product[END_REF]. This algorithm, intended to be used with rounding to nearest, is based on two error-free transformations: TwoSum [START_REF] Knuth | The Art of Computer Programming[END_REF] and TwoProd [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] that compute respectively the sum and the product of two floating-point numbers. The TwoProd algorithm requires 17 floating-point operations. However another error-free transformation, TwoProdFMA presented as Algorithm 7, exists for the product and costs only 2 floating-point operations ([18, p. 152]).

function [x, y] = TwoProdFMA(a, b) 1: x ← a × b 2: y ← FMA(a, b, -x) Algorithm 7:
Error-free transformation for the product of two floating-point numbers using an FMA

The TwoProdFMA algorithm is based on the Fused-Multiply-and-Add (FMA) operator that enables a floating-point multiplication followed by an addition to be performed as a single floating-point operation. For a, b, c ∈ F, FMA(a, b, c) is an approximation of a × b + c ∈ R that satisfies:

FMA(a, b, c) = (a × b + c)(1 + ε 1 ) = (a × b + c)/(1 + ε 2 )
where |ε ν | ≤ u with rounding to nearest and |ε ν | ≤ 2u with directed rounding. The FMA operation is supported by numerous processors such as AMD or Intel processors starting with respectively the Bulldozer or the Haswell architecture and by the Intel Xeon Phi coprocessor. It is also supported by AMD and NVidia GPUs (Graphics Processing Units) since 2010.

With any rounding mode, the TwoProdFMA algorithm computes both the floating-point product x of two numbers a and b and the associated rounding error y, provided that no underflow occurs. If this property holds, the floating-point numbers x and y computed by the TwoProdFMA algorithm satisfy:

x + y = a × b.
The CompDot algorithm, presented as Algorithm 8, is a compensated dot product algorithm based on FastTwoSum (Algorithm 2) and TwoProdFMA. As a remark, with rounding to nearest, the result of the CompDot algorithm is identical if other error-free transformations are used for the sum or the product. function res=CompDot(x, y)

1: [p, s] ← TwoProdFMA(x 1 , y 1 ) 2: for i = 2 to n do 3: [h, r] ← TwoProdFMA(x i , y i ) 4: [p, q] ← FastTwoSum(p, h) 5:
s ← s + (q + r) 6: end for 7: res ← p + s Algorithm 8: Compensated dot product of x = {x i } and y = {y i },

1 ≤ i ≤ n
The error on the result res of Algorithm 8 obtained with rounding to nearest is analysed in [START_REF] Ogita | Accurate sum and dot product[END_REF]. A bound for the absolute error is recalled in Proposition 5.3.

Proposition 5.3 ([20]

). Let floating-point numbers x i , y i ∈ F, 1 ≤ i ≤ n, be given and denote by res ∈ F the result computed by Algorithm 8 (CompDot) with rounding to nearest. If nu < 1, then,

|res -x T y| ≤ u|x T y| + γ 2 n (u)|x T ||y|.
(5.15)

In Corollary 5.4, Equation 5.15 is rewritten in terms of the condition number for the dot product.

Corollary 5.4 ([20]

). With rounding to nearest, if nu < 1, then, the result res of Algorithm 8

(CompDot) satisfies |res -x T y| |x T y| ≤ u + 1 2 γ 2 n (u) cond x T y .
As a consequence, we conclude that the result is as accurate as if computed in twice the working precision and then rounded to the current working precision. This is the same phenomenon as for the compensated summation algorithm.

Compensated dot product with directed rounding

We present here the impact of a directed rounding mode on Algorithm 8 (CompDot). For the error analysis we rewrite this algorithm into the following equivalent one. function res=CompDot(x, y)

1: [p 1 , s 1 ] ← TwoProdFMA(x 1 , y 1 ) 2: for i = 2 to n do 3: [h i , r i ] ← TwoProdFMA(x i , y i ) 4: [p i , q i ] ← FastTwoSum(p i-1 , h i ) 5: s i ← s i-1 + (q i + r i ) 6: end for 7: res ← p n + s n Algorithm 9: Equivalent formulation of Algorithm 8
A bound for the absolute error on the result res of Algorithm 8 is given in Proposition 5.5.

Proposition 5.5. Let floating-point numbers x i , y i ∈ F, 1 ≤ i ≤ n, be given and denote by res ∈ F the result computed by Algorithm 8 (CompDot) with directed rounding. If (n + 1)u < 1 2 , then,

|res -x T y| ≤ 2u|x T y| + 2γ 2 n+1 (2u)|x T ||y|.
Proof. Thanks to the TwoProdFMA algorithm, we have

p 1 + s 1 = x 1 y 1 , (5.16)
and for i ≥ 2, h i + r i = x i y i .

(5.17)

From Proposition 4.5, it follows that

p i + e i = p i-1 + h i with |q i -e i | ≤ 2u|e i |.
(5.18)

Therefore from Equation 5.17, we deduce that

e i + r i = (p i-1 + h i -p i ) + (x i y i -h i ) = x i y i + p i-1 -p i .
Then from Equation 5.16, we derive

s 1 + n i=2 (e i + r i ) = (x 1 y 1 -p 1 ) + n i=2 x i y i + p 1 -p n = x T y -p n . (5.19)
Because the TwoProdFMA algorithm is executed with a directed rounding mode, for i ≥ 2, then

|r i | ≤ 2u|x i y i |.
Therefore, we have

n i=2 |r i | ≤ 2u n i=2 |x i y i |,
and

|s 1 | + n i=2 |r i | ≤ 2u|x T ||y|.
(5.20)

From Lemma 4.6, we deduce

n i=2 |e i | ≤ γ n-1 (2u) |p 1 | + n i=2 |h i | .
As a consequence, we have

n i=2 |e i | ≤ γ n-1 (2u) n i=1 | fl * (x i y i )| ,
and

n i=2 |e i | ≤ (1 + 2u)γ n-1 (2u)|x T ||y|. ( 5
|q i | ≤ n i=2 |e i | + n i=2 |q i -e i | .
Therefore, Equation 5.23 yields

n i=2 |q i | ≤ (1 + 2u) n i=2 |e i | .
From Equations 4.12 and 5.22, it yields

n i=2 |q i | ≤ γ n+1 (2u)|x T ||y|. (5.25)
For later use, we evaluate an upper bound on the following expression

|s 1 + n i=2 (q i + r i ) -s n | = s 1 + n i=2 (q i + r i ) -fl * s 1 + n i=2 (q i + r i ) .
From Proposition 4.1, it follows that

|s 1 + n i=2 (q i + r i ) -s n | ≤ γ n-1 (2u) |s 1 | + n i=2 |fl * (q i + r i )| . (5.26)
Furthermore, because a directed rounding mode is used, we have

n i=2 |fl * (q i + r i )| ≤ (1 + 2u) n i=2 |q i + r i | .
Therefore from Equation 5.26, we deduce that

|s 1 + n i=2 (q i + r i ) -s n | ≤ (1 + 2u)γ n-1 (2u) |s 1 | + n i=2 |q i + r i | ,
and, from Equation 4.12,

|s 1 + n i=2 (q i + r i ) -s n | ≤ γ n (2u) |s 1 | + n i=2 |q i + r i | .
From Equations 5.20 and 5.25, it follows that

|s 1 + n i=2 (q i + r i ) -s n | ≤ γ n (2u) (2u + γ n+1 (2u)) |x T ||y|.
(5.27)

We deduce from Equation 5.19 that

| x T y -p n -s n | = s 1 + n i=2 (e i + r i ) -s n .
As a consequence, it yields

|x T y -p n -s n | = s 1 + n i=2 (q i + r i ) -s n + n i=2 (e i -q i ) ,
and

|x T y -p n -s n | ≤ s 1 + n i=2 (q i + r i ) -s n + n i=2 |e i -q i | .
Therefore from Equations 5.24 and 5.27, we deduce that

|x T y -p n -s n | ≤ γ n (2u) (4u + γ n+1 (2u)) |x T ||y|.
(5.28)

Let us show that γ n+1 (2u) ≥ 4u. It is easy to show that

γ n+1 (2u) -4u = 2(n + 1)u 1 -2(n + 1)u -4u,
and γ n+1 (2u) -4u = 2u 1 -2(n + 1)u (n -1 + 4(n + 1)u) .

Because (n + 1)u < 1 2 , it follows that γ n+1 (2u) -4u ≥ 0.

Therefore from Equation 5.28, we can deduce that

|x T y -p n -s n | ≤ 2γ n (2u)γ n+1 (2u)|x T ||y|. (5.29)
Because Algorithm 9 is executed with a directed rounding mode, it follows that

|res -x T y| = |(1 + ε)(p n + s n ) -x T y| with |ε| ≤ 2u.
Therefore, we have

|res -x T y| = |εx T y + (1 + ε)(p n + s n -x T y)|, and |res -x T y| ≤ 2u|x T y| + (1 + 2u)|p n + s n -x T y|.
Then from Equation 5.29, it follows that

|res -x T y| ≤ 2u|x T y| + 2(1 + 2u)γ n (2u)γ n+1 (2u)|x T ||y|.
Finally from Equation 4.12, we conclude that

|res -x T y| ≤ 2u|x T y| + 2γ 2 n+1 (2u)|x T ||y|.
From Proposition 5.5, a bound for the relative error on the result of Algorithm 8 (CompDot) obtained with directed rounding is deduced in Corollary 5.6.

Corollary 5.6. With directed rounding, if (n + 1)u < 1 2 , then, the result res of Algorithm 8 (CompDot) satisfies |res -

x T y| |x T y| ≤ 2u + γ 2 n+1 (2u) cond x T y .
From Corollary 5.6, the relative error bound on the result of Algorithm 8 (CompDot) computed with directed rounding is essentially (nu) 2 times the condition number plus the rounding 2u due to the working precision. Like with rounding to nearest, the result obtained with directed rounding is as accurate as if computed in twice the working precision and then rounded to the current working precision.

Accurate Horner scheme

In this section, we present the accuracy obtained with the classic Horner scheme for polynomial evaluation. We recall a compensated Horner scheme algorithm and the error on its result computed with rounding to nearest. Then we analyse the impact of a directed rounding mode on this algorithm. In this section, we assume that no underflow occurs.

Classic Horner scheme

The classical method for evaluating a polynomial If an FMA instruction is available, then the statement s i ← s i+1 • x + a i in Algorithm 10 can be rewritten as s i ← FMA(s i+1 , x, a i ) which slightly improves the error bound. Using an FMA this way, the computed result now satisfies |p(x) -res| ≤ γ n (2u) p(|x|).

A compensated Horner scheme with rounding to nearest

We now want to accurately compute a polynomial at a given point. The Horner scheme algorithm can be modified to compute the rounding error at each elementary operation using error-free transformations. We present as Algorithm 11 a compensated algorithm for the Horner scheme. One can find a more detailed description of the compensated Horner scheme algorithm in [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF][START_REF] Graillat | Compensated Horner scheme[END_REF]. Algorithm 11 is based on the TwoProdFMA and FastTwoSum algorithms. However, with rounding to nearest, its result is identical if other error-free transformations, such as TwoProd or TwoSum, are used. function res = CompHorner(p, x)

1: s n ← a n 2: r n ← 0 3: for i = n -1 down to 0 do 4: [p i , π i ] ← TwoProdFMA(s i+1 , x) 5: [s i , σ i ] ← FastTwoSum(p i , a i ) 6: r i ← r i+1 • x + (π i + σ i ) 7:
end for 8: res ← s 0 + r 0 Algorithm 11: Polynomial evaluation with a compensated Horner scheme If we denote by p π and p σ the two following polynomials

p π (x) = n-1 i=0 π i x i , p σ (x) = n-1 i=0 σ i x i ,
then one can show, thanks to error-free transformations, that

p(x) = s 0 + p π (x) + p σ (x).
If one looks closely at the previous algorithm, it is then clear that s 0 = Horner(p, x). As a consequence, we can derive a new error-free transformation for the polynomial evaluation p(x) = Horner(p, x) + p π (x) + p σ (x).

The compensated Horner scheme first computes p π (x) + p σ (x) which corresponds to the rounding errors and then adds the value obtained to the result of the classic Horner scheme Horner(p, x). We will show that the result computed by Algorithm 11 admits a significantly better error bound than that computed with the classical Horner scheme. We argue that Algorithm 11 provides a result as if it was computed using twice the working precision. This is summed up in the following theorem. In other words, the bound for the relative error of the computed result is essentially γ 2 2n (u) times the condition number of the polynomial evaluation, plus the unavoidable term u for rounding the result to the working precision. In particular, if cond(p, x) < γ -1 2n (u), then the relative accuracy of the result is bounded by a constant of the order of u. This means that the compensated Horner algorithm computes an evaluation accurate to the last few bits as long as the condition number is smaller than γ -1 2n (u) ≈ (2nu) -1 . Besides that, (6.32) tells us that the computed result is as accurate as if computed by the classic Horner algorithm with twice the working precision, and then rounded to the working precision.

A compensated Horner scheme with directed rounding

We now present the impact of a directed rounding mode on Algorithm 11 (CompHorner).

Let τ i be the rounding error in the floating-point addition of p i and a i (τ i is not necessarily a floating-point number):

s i + τ i = p i + a i .
It follows that s i+1 • x = p i + π i and p i + a i = s i + τ i with |τ i -σ i | ≤ 2uτ i . As a consequence, we have

s i = s i+1 • x -π i -τ i for i = 0, . . . , n -1.
By induction, we deduce that p(x) = s 0 + p π (x) + p τ (x), with

s 0 = fl * (p(x)), p π (x) = n-1 i=0 π i x i , and p τ (x) = n-1 i=0 τ i x i . (6.33)
We recall that

p σ (x) = n-1 i=0 σ i x i . (6.34)
In the sequel, we will denote e(x) := p π (x) + p σ (x). In this case, we have p(x) = fl(p(x)) + e(x) + (p τ -p σ )(x) and res = fl(p(x) + e(x)). Lemma 6.2. Let p(x) = n i=0 a i x i a polynomial with a i ∈ F, 0 ≤ i ≤ n and x ∈ F. Let p π and p σ be defined by (6.33) and (6.34). Then, we have

p π (|x|) + p σ (|x|) ≤ γ 2n+1 (2u) p(|x|), with p(x) = n i=0 |a i |x i .
Proof. Using (2.1), we have, for i = 1, . . . , n,

|p n-i | = | fl * (s n-i+1 • x)| ≤ (1 + 2u)|s n-i+1 ||x| and |s n-i | = | fl * (p n-i + a n-i )| ≤ (1 + 2u)(|p n-i | + |a n-i |).
Let us show by induction on i = 1, . . . , n that

|p n-i | ≤ (1 + γ 2i-1 (2u)) i j=1 |a n-i+j ||x j | (6.35)
and

|s n-i | ≤ (1 + γ 2i (2u)) i j=0 |a n-i+j ||x j |. (6.36)
For i = 1, as s n = a n , we have

|p n-1 | ≤ (1 + 2u)|a n ||x| ≤ (1 + γ 1 (2u))|a n ||x|
and so (6.35) is true. In the same way, as

|s n-1 | ≤ (1 + 2u) ((1 + γ 1 (2u))|a n ||x| + |a n-1 |) ≤ (1 + γ 2 (2u)) (|a n ||x| + |a n-1 |)
then (6.36) is also true. Let us assume that (6.35) and (6.36) are true for an integer i, 1 ≤ i < n.

Then we have

|p n-(i+1) | ≤ (1 + 2u)|s n-i ||x|.
By hypothesis, we deduce that

|p n-(i+1) | ≤ (1 + 2u)(1 + γ 2i (2u)) i j=0 |a n-i+j ||x j+1 | ≤ (1 + γ 2(i+1)-1 (2u)) i+1 j=1 |a n-(i+1)+j ||x j |.
Hence, it follows that

|s n-(i+1) | ≤ (1 + 2u)(|p n-(i+1) | + |a n-(i+1) |) ≤ (1 + 2u)(1 + γ 2(i+1)-1 (2u))   i+1 j=1 |a n-(i+1)+j ||x j | + |a n-(i+1) |   ≤ (1 + γ 2(i+1) (2u)) i+1 j=0 |a n-(i+1)+j ||x j |.
Relations (6.35) and (6.36) are then true by induction. As a consequence, for i = 1, . . . , n, we have

|p n-i ||x n-i | ≤ (1 + γ 2i-1 (2u)) p(x)
and

|s n-i ||x n-i | ≤ (1 + γ 2i (2u)) p(x).

Following (2.2), we have |π

i | ≤ 2u|p i |, |τ i | ≤ 2u|s i | and |σ i | ≤ (1 + 2u)|τ i | for i = 0, . . . , n -1. Hence, ( p π + p σ )(|x|) = n-1 i=0 (|π i | + |σ i |)|x i | ≤ 2u(1 + 2u) n i=1 (|p n-i | + |s n-i |)|x n-i |,
and so

( p π + p σ )(|x|) ≤ 2u(1 + 2u) n i=1 (2 + γ 2i-1 (2u) + γ 2i (2u)) p(|x|) ≤ 4nu(1 + 2u) (1 + γ 2n (2u)) p(|x|). As 4nu(1 + γ 2n (2u)) = γ 2n (2u), we deduce that ( p π + p σ )(|x|) ≤ γ 2n+1 (2u) p(|x|). Lemma 6.3. Let p(x) = n i=0 a i x i be a polynomial with a i ∈ F, 0 ≤ i ≤ n, q(x) = n i=0 b i x i a polynomial with b i ∈ F, 0 ≤ i ≤ n and x ∈ F.
Then the floating-point evaluation of r(x) = p(x) + q(x) via the following algorithm

1: r n ← fl * (a n + b n ) 2: for i = n -1 down to 0 do 3: r i ← fl * (r i+1 • x + (a i + b i )) 4: end for 5: res ← r 0 satisfies |res -r(x)| ≤ γ 2n+1 (2u) r(|x|).
Proof. Considering the previous algorithm, we have r n = fl * (a n + b n ) = (a n + b n ) 1 (2u) and for i = n -1 down to 0,

r i = fl * (r i+1 • x + (a i + b i )) = r i+1 x 2 (2u) + (a i + b i ) 2 (2u).
As a consequence, we can show by induction that

r 0 = (a n + b n )x n 2n + 1 (2u) + n-1 i=0 (a i + b i )x i 2(i + 1) (2u).
Moreover with quantities θ 2n+1 (2u), θ 2n (2u), . . . , θ 1 (2u), satisfying |θ i (2u)| ≤ γ i (2u), we have

r 0 = (a n + b n )x n (1 + θ 2n+1 )(2u) + n-1 i=0 (a i + b i )x i (1 + θ 2(i+1) (2u)).
As r 0 = fl * (p(x) + q(x)), we finally get

res - n i=0 (a i + b i )x i ≤ γ 2n+1 (2u) n i=0 |a i + b i ||x i | ≤ γ 2n+1 (2u)( p + q)(|x|),
which concludes the proof. Proof. By considering Algorithm 11, we have p(x) = s 0 + e(x) + (p τ -p σ )(x). We can deduce that

|res -p(x)| = |(1 + ε)(s 0 + fl * (e(x))) -p(x)| = |(1 + ε)(s 0 + fl * (e(x)) -p(x) + (p τ -p σ )(x)) + εp(x) + (1 + ε)(p σ -p τ )(x)| = |(1 + ε)(s 0 + e(x) + (p τ -p σ )(x) -p(x)) + (1 + ε)(fl * (e(x)) -e(x))+ εp(x) + (1 + ε)(p τ -p σ )(x)| ≤ 2u|p(x)| + (1 + 2u)| fl * (e(x)) -e(x)| + (1 + 2u)|(p τ -p σ )(x)|.
By applying Lemma 6.3, we obtain

| fl * (e(x)) -e(x)| ≤ γ 2n-1 (2u) e(|x|) ≤ γ 2n-1 (2u)( p π (|x|) + p σ (|x|)).
Moreover from Lemma 6.2, we get As (1 + 2u)γ 2n-1 (2u) ≤ γ 2n (2u) and 2nu ≤ γ 2n+1 (2u), we obtain

p π (|x|) + p σ (|x|) ≤ γ 2n+1 (2u) p(|x|). Since |τ i -σ i | ≤ 2uτ i , we have |(p τ -p σ )(x)| ≤ 2u n-1 i=0 |τ i ||x| i ≤ 2u p τ (|x|).
|res -p(x)| ≤ 2u|p(x)| + 2γ 2n+1 (2u) 2 p(|x|),
which concludes the proof.

Summation as in K-fold precision

According to Sect. 4, Algorithms 3 (FastCompSum) and 5 (PriestCompSum) compute the sum of n floating-point numbers as in twice the working precision, even with directed rounding. In this section, we present the SumK algorithm [START_REF] Ogita | Accurate sum and dot product[END_REF] that computes this sum as in K-fold precision. We recall the error bound on its result obtained with rounding to nearest. Then we analyse the impact of directed rounding on this algorithm.

Summation as in K-fold precision with rounding to nearest

The SumK algorithm, introduced in [START_REF] Ogita | Accurate sum and dot product[END_REF] and presented as Algorithm 12, enables the summation of a vector of floating-point numbers as in K-fold precision. In [START_REF] Ogita | Accurate sum and dot product[END_REF] the SumK algorithm is based on the TwoSum algorithm [START_REF] Knuth | The Art of Computer Programming[END_REF]. However the same result res can be obtained using rounding to nearest with any error-free transformation that computes the sum of two floating-point numbers. The SumK algorithm is presented here with Algorithm 4 (PriestTwoSum) that, although more costly, is an error-free transformation with any rounding mode. As a remark if K = 2, Algorithm 12 is identical with Algorithm 5 (PriestCompSum).

function res = SumK(p, K)

1: for k = 1 to K -1 do 2:
for i = 2 to n do 3:

[p i , p i-1 ] ← PriestTwoSum(p i , p i-1 ) 4:
end for 5: end for 6: res ← n i=1 p i Algorithm 12: Summation of n floating-point numbers p = {p i } in K-fold working precision, K ≥ 3

The error on the result res of Algorithm 12 obtained with rounding to nearest is analysed in [START_REF] Ogita | Accurate sum and dot product[END_REF]. A bound for the absolute error is recalled in Proposition 7.1 and a bound for the relative error in Corollary 7.2.

Proposition 7.1 ([20]

). Let floating-point numbers p i ∈ F, 1 ≤ i ≤ n, be given and assume 4nu ≤ 1. Then, also in the presence of underflow, the result res of Algorithm 12 (SumK) obtained with rounding to nearest satisfies for K ≥ 3

|res -s| ≤ u + 3γ 2 n-1 (u) |s| + γ K 2n-2 (u)S
where s := p i and S := |p i |.

Corollary 7.2 ([20]

). Assume 4nu ≤ 1. The result res of Algorithm 12 (SumK) obtained with rounding to nearest, also in the presence of underflow, satisfies

|res -s| |s| ≤ u + 3γ 2 n-1 (u) + γ K 2n-2 (u) cond p i .
From Corollary 7.2, because γ n (u) ≈ nu, the bound for the relative error on the result obtained with rounding to nearest is essentially the relative rounding error u plus a term that reflects that the computation is carried out as in K-fold precision ((αu) K times the condition number for a moderate factor α).

Summation as in K-fold precision with directed rounding

We analyse here the impact of a directed rounding mode on Algorithm 12 (SumK). The steps of the error analysis are similar to those presented in [START_REF] Ogita | Accurate sum and dot product[END_REF] for rounding to nearest. We denote the input vector p by p (0) , and the vector after finishing loop k by p (k) . We also set

S (k) := n i=1 |p (k) i | for 0 ≤ k ≤ K -1.
s := n i=1 p (0) i = n i=1 p (k) i for 1 ≤ k ≤ K -1, (7.37) |res -s| ≤ 2u|s| + γ 2 n-1 (2u)S (K-2) , (7.38)
S (k) ≤ 3|s| + γ k 2n-2 (2u)S (0) provided 8(n -1)u ≤ 1 and 1 ≤ k ≤ K -1. (7.39)
Proof. Equation 7.37 follows by successive applications of Equation 4.10. Equation 7.38 can be deduced using s = n i=1 p (K-2) i and applying Proposition 4.10. Let us now prove Equation 7.39. From Lemma 4.9, we obtain

n i=1 |p (1) i | ≤ |s| + 2γ n-1 (2u)S (0) . (7.40)
By applying successively Equation 7.40 and using Equation 7.37 we obtain

S (2) ≤ |s| + 2γ n-1 (2u) |s| + 2γ n-1 (2u)S (0) , and
S (k) ≤ |s| + ∞ i=0 (2γ n-1 (2u)) i + (2γ n-1 (2u)) k S (0) for 1 ≤ k ≤ K -1.
We have

∞ i=0 (2γ n-1 (2u)) i = 1 -2(n -1)u 1 -6(n -1)u . If 8(n -1)u ≤ 1, then we get 1 -2(n -1)u ≤ 3(1 -6(n -1)u) and 1 -2(n -1)u 1 -6(n -1)u ≤ 3.
Therefore, we have

S (k) ≤ 3|s| + (2γ n-1 (2u)) k S (0) .
Because 2γ m (2u) ≤ γ 2m (2u), we can conclude that

S (k) ≤ 3|s| + (γ 2n-2 (2u)) k S (0) .
A bound for the absolute error on the result of Algorithm 12 (SumK) obtained with directed rounding is given in Proposition 7.4. Proposition 7.4. Let floating-point numbers p i ∈ F, 1 ≤ i ≤ n, be given and assume 8nu ≤ 1. Then, also in the presence of underflow, the result res of Algorithm 12 (SumK) obtained with directed rounding satisfies for K ≥ 3

|res -s| ≤ 2u + 3γ 2 n-1 (2u) |s| + γ K 2n-2 (2u)S
where s := p i and S := |p i |.

Proof. Proposition 7.4 can be proved by inserting Equation 7.39 of Lemma 7.3 into Equation 7.38.

From Proposition 7.4, a bound for the relative error on the result of Algorithm 12 (SumK) obtained with directed rounding is deduced in Corollary 7.5.

Corollary 7.5. Assume 8nu ≤ 1. The result res of Algorithm 12 (SumK) obtained with directed rounding, also in the presence of underflow, satisfies

|res -s| |s| ≤ 2u + 3γ 2 n-1 (2u) + γ K 2n-2 (2u) cond p i . (7.41)
From Corollary 7.5, because γ n (2u) ≈ 2nu, the bound for the relative error on the result res obtained with directed rounding is essentially the relative rounding error 2u plus (αu) K times the condition number for a moderate factor α. Like with rounding to nearest, the last term on the right hand side of Equation 7.41 reflects that the computation is carried out as in K-fold precision.

8 Dot product as in K-fold precision According to Sect. 5, Algorithm 8 (CompDot) computes a dot product as in twice the working precision, even with directed rounding. In this section, we present the DotK algorithm [START_REF] Ogita | Accurate sum and dot product[END_REF] that computes a dot product as in K-fold precision. We recall the error bound on its result obtained with rounding to nearest. Then we analyse the impact of directed rounding on this algorithm. Like in Sect. 5, we assume in this section that no underflow occurs.

Dot product as in K-fold precision with rounding to nearest

The DotK algorithm, introduced in [START_REF] Ogita | Accurate sum and dot product[END_REF] and presented as Algorithm 13, enables one to compute a dot product as in K-fold precision. In [START_REF] Ogita | Accurate sum and dot product[END_REF] the DotK algorithm is based on TwoProd [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] and TwoSum [START_REF] Knuth | The Art of Computer Programming[END_REF] that are error-free transformations with rounding to nearest. The DotK algorithm is presented here with TwoProdFMA and PriestTwoSum that are error-free transformations with any rounding mode. As a remark if K = 2, Algorithm 13 is identical with Algorithm 8 (CompDot). function res = DotK(x, y, K)

1: [p, r 1 ] ← TwoProdFMA(x 1 , y 1 ) 2: for i = 2 to n do 3: [h, r i ] ← TwoProdFMA(x i , y i ) 4:
[p, r n+i-1 ] ← PriestTwoSum(p, h) 5: end for 6: r 2n ← p 7: res ← SumK(r, K -1)

Algorithm 13: Dot product algorithm in K-fold working precision, K ≥ 3

The error on the result res of Algorithm 13 obtained with rounding to nearest is analysed in [START_REF] Ogita | Accurate sum and dot product[END_REF]. A bound for the absolute error is recalled in Proposition 8. [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF]

-x T y x T y ≤ u + 2γ 2 4n-2 (u) + 1 2 γ K 4n-2 (u) cond x T y .
From Corollary 8.2, the bound for the relative error on the result is essentially the relative rounding error u plus a term that reflects that the computation is carried out as in K-fold precision (α(K)u K times the condition number for a moderate factor α(K)).

Dot product as in K-fold precision with directed rounding

We present here the impact of a directed rounding mode on Algorithm 13 (DotK). For the analysis we rewrite Algorithm 13 into the following equivalent one. function res = DotK(x, y, K) Indeed, it is clear that

1: [p 1 , r 1 ] ← TwoProdFMA(x 1 , y 1 ) 2: for i = 2 to n do 3: [h i , r i ] ← TwoProdFMA(x i , y i ) 4: [p i , r n+i-1 ] ← PriestTwoSum(p i-1 , h i
r 1 = x 1 y 1 -p 1 ,
and for i ≥ 2,

r i + r n+i-1 = (x i y i -h i ) + (p i-1 + h i -p i ), = x i y i + p i-1 -p i .
Therefore, because r 2n = p n , we have

2n-1 i=1 r i = (x 1 y 1 -p 1 ) + n i=2 x i y i + p 1 -r 2n , = x T y -r 2n . ( 8
|r n+i-1 | ≤ γ n-1 (2u) |p 1 | + n i=2 |h i | , = γ n-1 (2u) n i=1 |f l * (x i y i )|, ≤ (1 + 2u)γ n-1 (2u)|x T ||y|. ( 8
|r i | ≤ 2u|x T ||y| + (1 + 2u)γ n-1 (2u)|x T ||y|, ≤ 2nu 1 -2(n -1)u |x T ||y|. (8.47)
From Equation 8.43, we deduce that

|r 2n | = |x T y - 2n-1 i=1 r i |, ≤ |x T y| + 2n-1 i=1 |r i |.
Therefore, we have

2n i=1 |r i | ≤ |x T y| + 2 2n-1 i=1 |r i |. (8.48)
From Equation 8.47, we obtain 2 2n-1 From Proposition 7.4, using 2γ m (2u) ≤ γ 2m (2u) and noting that the vector r is of length 2n yields

i=1 |r i | ≤ (2n)(2u) 1 -(n -1)(2u) |x T ||y|, ≤ (2n)
|res -x T y| ≤ 2u + 3γ 2 2n-1 (2u) |x T y| + γ K-1 4n-2 (2u) |x T y| + γ 2n (2u)|x T ||y| , ≤ 2u + 3γ 2 2n-1 (2u) + γ K-1 4n-2 (2u) |x T y| + γ 2n (2u)γ K-1 4n-2 (2u)|x T ||y|, ≤ 2u + 3 4 γ 2 4n-2 (2u) + γ K-1 4n-2 (2u) |x T y| + γ K 4n-2 (2u)|x T ||y|.
If 8(2n -1)u ≤ 1, then we can conclude that γ 4n-2 (2u) ≤ 1 and From Corollary 8.4, the bound for the relative error on the result obtained with directed rounding is essentially the relative rounding error 2u plus a term that reflects that the computation is carried out as in K-fold precision (α(K)u K times the condition number for a moderate factor α(K)).

|res -x T y| ≤ 2u + 2γ 2 4n-2 (2u) |x T y| + γ K 4n-2 (

Numerical results

In the numerical experiments presented here, compensated algorithms previously described are executed with the CADNA library. Results are computed in double precision, i.e. using the binary64 format of the IEEE standard [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF]: each stochastic variable contains three binary64 floatingpoint values. Figures 1 to 5 present the number of exact significant digits estimated by CADNA in the results obtained using classic algorithms and associated compensated versions. One can observe in Fig. 1 to 5 that, if the condition number increases, the number of exact significant digits decreases and with classic algorithms, results have no more correct digits for condition numbers greater than 10 16 .

From Fig. 1 to 3, as long as the condition number is less than 10 16 , the compensated algorithms produce results with the maximal accuracy (15 exact significant digits in double precision). For condition numbers greater than 10 16 , the accuracy decreases and there are no more correct digits for condition numbers greater than 10 32 . The results provided by CADNA are consistent with the properties of compensated algorithms given in Sect. 4 to 6 for directed rounding: with the current precision, the FastCompSum, CompDot, and CompHorner algorithms compute results that could have been obtained with twice the working precision. the accuracy estimated by CADNA of the results computed using the SumK and DotK algorithms described in Sect. 7 and 8. As a remark, although the results obtained with SumK and DotK for K = 2 are reported in Fig. 4 and5, for performance reasons algorithms FastCompSum and CompDot should be preferably used for a computation as with twice the working precision. It can be observed in Fig. 4 and 5 that, if the condition number is less than about 10 16(K-1) , algorithms SumK and DotK produce results with the maximum possible accuracy in double precision. Then, the accuracy decreases if the condition number increases from about 10 16(K-1) to 10 16K . For condition numbers greater than about 10 16K , results computed by SumK and DotK have no more correct digits. Figures 4 and5 CADNA can detect numerical instabilities that occur during the execution. As already mentioned in Sect. 3, the use of DSA requires self-validation, i.e. the control of multiplications and divisions duing the execution. Indeed if each operand of a multiplication or a divisor has no more correct digits, the estimation of accuracy by DSA may be invalid. The algorithms presented in this article require no division. No multiplication is performed in summation algorithms. Concerning dot product algorithms, the operands of multiplications are elements of the initial arrays. For the evaluation of a polynomial p(x), the value x is always one of the operands of the multiplications performed. Because these initial data are supposed to be significant (not numerical noise), no unstable multiplication can be generated. All algorithms presented in this article may generate cancellations, i.e. sudden losses of accuracy due to subtractions of very close values. The number of cancellations detected depends on the condition number, the size of the data arrays for the sum and the dot product, and the degree of the polynomial for Horner scheme. Because cancellations are not related to the self-validation of DSA, they cannot invalid the estimation of accuracy by CADNA.

Tables 1 to 3 present execution times measured in double precision with and without CADNA on an Intel Core i5-4690 CPU (Haswell) at 3.5 GHz using g++ version 5.3.1. As a remark, this architecture supports the FMA operation. The codes have been run using CADNA with three kinds of instability detection:

• no detection of instabilities,

• self-validation,

• the detection of all kinds of instabilities;

With the algorithms considered in this article, the execution times measured with self-validation are very close to those obtained if instability detection is deactivated. With summation algorithms, these times are necessarily the same. Therefore the execution times reported in Tables 1 to 3 have been measured with self-validation or with the detection of all kinds of instabilities.

From Tables 1 to 3, the cost of compensated algorithms that compute results as with twice the working precision (FastCompSum, CompDot, CompHorner) over the classic algorithms is about 2 without CADNA and about 3 with CADNA. The heavier cost of compensated algorithms with CADNA is mainly explained by the increase of data movements that are more costly with stochastic variables. The executions times of PriestCompSum and SumK for K = 2 are mentioned in Table 1. However, for performance reasons, FastCompSum that is based on FastTwoSum should be preferably used to compute summations as with twice the working precision. Indeed PriestCompSum is based on PriestTwoSum that requires more floating-point operations and an extra branching statement compared to FastTwoSum. SumK for K = 2 is more costly than PriestCompSum because it requires to fetch the errors from an array. Indeed at the end of the computation SumK sums up the errors stored in an array, while PriestCompSum adds each error as soon as it is available. Similarly the execution times of DotK for K = 2 are reported in Table 2, although CompDot should be preferable used to compute dot products as with twice the working precision. Indeed PriestTwoSum that is called in DotK is more costly than FastTwoSum called in CompDot. Furthermore DotK for K = 2 sums up the errors fetched from an array at the end of the computation, while CompDot adds the errors as soon as they are computed. One can observe in Table 1 that the cost of SumK over the classic summation regularly increases with K, whatever the instability detection level with CADNA and also without CADNA. A similar remark can be formulated from Table 2 about the cost of DotK over the classic dot product. In all the algorithms mentioned in Tables 1 to 3, the cost of CADNA

Conclusion

We have shown that it is possible to validate some compensated algorithms using stochastic arithmetic. We studied compensated summation, dot product and polynomial evaluation. For that, we described the behavior of error-free transformations for addition and multiplication when a random rounding mode is used. We were also able to validate K-fold compensated algorithms for summation and dot-product by using a special error-free transformation from Priest [START_REF] Priest | On Properties of Floating Point Arithmetics: Numerical Stability and the Cost of Accurate Computations[END_REF].

As a future work, we plan to show that it may be possible to validate other compensated algorithms with stochastic arithmetic. We intend to study compensated floating-point product and exponentiation [START_REF] Graillat | Accurate floating point product and exponentiation[END_REF], compensated Newton's scheme [START_REF] Graillat | Accurate simple zeros of polynomials in floating point arithmetic[END_REF][START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial[END_REF], and the compensated evaluation of elementary symmetric functions [START_REF] Jiang | Accurate and fast evaluation of elementary symmetric functions[END_REF]. Moreover, a more challenging problem will be to validate the compensated algorithm for solving triangular systems [START_REF] Louvet | Algorithmes compensés en arithmétique flottante : précision, validation, performances[END_REF].

. 1 )

 1 As a consequence,|a • b -fl * (a • b)| ≤ 2u|a • b| and |a • b -fl * (a • b)| ≤ 2u| fl * (a • b)| for • = {+, -}.(2.2)

  3 and a bound for the relative error in Corollary 4.4. function [c, d] = FastTwoSum(a, b) 1: if |b| > |a| then 2: exchange a and b 3: end if 4: c ← a + b 5: z ← c -a 6: d ← b -z Algorithm 2: Error-free transformation for the sum of two floating-point numbers with rounding to nearest function res = FastCompSum(p)

Proposition 4 . 5 (

 45 [START_REF] Graillat | Numerical validation of compensated summation algorithms with stochastic arithmetic[END_REF]). Let c and d be the floating-point addition of a and b and the correction both computed by Algorithm 2 (FastTwoSum) using directed rounding. Let e be the error on c: a + b = c + e. Then |e -d| ≤ 2u|e|.

  then

4 :

 4 Error-free transformation for the sum of two floating-point numbers with any rounding mode A compensated summation algorithm based on PriestTwoSum is given as Algorithm 5.

  i is the Horner scheme which consists of Algorithm 10. function res = Horner(p, x) 1: s n ← a n 2: for i = n -1 downto 0 do 3: s i ← s i+1 • x + a i 4: end for 5: res ← s 0 Algorithm 10: Polynomial evaluation with Horner's scheme Whatever the rounding mode, a forward error bound on the result of Algorithm 10 is (see [11, p. 95]): |p(x) -res| ≤ γ 2n (2u) n i=0 |a i ||x| i = γ 2n (2u) p(|x|) where p(x) = n i=0 |a i |x i . It is very interesting to express and interpret this result in terms of the condition number of the polynomial evaluation defined by cond(p, x) = n i=0 |a i ||x| i |p(x)| = p(|x|) |p(x)| . (6.30) Thus we have |p(x) -res| |p(x)| ≤ γ 2n (2u) cond(p, x).

Theorem 6 . 4 .

 64 Consider a polynomial p of degree n with floating-point coefficients, and a floatingpoint value x. With directed rounding, the forward error in the compensated Horner algorithm is such that |CompHorner(p, x) -p(x)| ≤ 2u|p(x)| + 2γ 2n+1 (2u) 2 p(x).

Moreover, as |τ

  i | ≤ 2u|s i |, we have p τ (|x|) ≤ 2nuγ 2n (2u) p(|x|). As a consequence, we deduce |res -p(x)| ≤ 2u|p(x)| + (1 + 2u)γ 2n-1 (2u)γ 2n+1 (2u) p(|x|) + 2nu(1 + 2u)γ 2n (2u) p(|x|).

Lemma 7 . 3 .

 73 With the above notations, the intermediate results of Algorithm 12 (SumK) satisfy:

) 5 :

 5 end for 6: r 2n ← p n 7: res ← SumK(r, K -1) Algorithm 14: Equivalent formulation of Algorithm 13 A bound for the absolute error on the result of Algorithm DotK obtained with directed rounding is given in Proposition 8.3. Proposition 8.3. Let floating-point numbers x i , y i ∈ F, 1 ≤ i ≤ n, be given and denote by res ∈ F the result computed by Algorithm 13 (DotK) with directed rounding. If 16nu ≤ 1, then |res -x T y| ≤ 2u + 2γ 2 4n-2 (2u) |x T y| + γ K 4n-2 (2u)|x T ||y| Proof. Because TwoProdFMA and PriestTwoSum are error-free transformations even with directed rounding, we have s := 2n i=1 r i = x T y. (8.42)

Figure 1 :

 1 Figure 1: Accuracy estimated by CADNA using the Sum and the FastCompSum algorithms with 200 randomly generated floating-point numbers.

Figure 2 :

 2 Figure 2: Accuracy estimated by CADNA using the Dot and the CompDot algorithms to compute the dot product of arrays of 100 randomly generated elements.

Figure 3 :

 3 Figure 3: Accuracy estimated by CADNA using the Horner and the CompHorner algorithms to compute (x -1) n for x close to 1 and for various values of n.

Figures 4 and 5

 5 Figures 4 and 5 present the accuracy estimated by CADNA of the results computed using the SumK and DotK algorithms described in Sect. 7 and 8. As a remark, although the results obtained with SumK and DotK for K = 2 are reported in Fig.4 and 5, for performance reasons algorithms FastCompSum and CompDot should be preferably used for a computation as with twice the working precision. It can be observed in Fig.4and 5 that, if the condition number is less than about 10 16(K-1) , algorithms SumK and DotK produce results with the maximum possible accuracy in double precision. Then, the accuracy decreases if the condition number increases from about 10 16(K-1) to 10 16K . For condition numbers greater than about 10 16K , results computed by SumK and DotK have no more correct digits. Figures4 and 5are consistent with the properties of algorithms SumK and DotK given in Sect. 7 and 8 for directed rounding: results are computed as in K-fold precision.

  Figures 4 and 5 present the accuracy estimated by CADNA of the results computed using the SumK and DotK algorithms described in Sect. 7 and 8. As a remark, although the results obtained with SumK and DotK for K = 2 are reported in Fig.4 and 5, for performance reasons algorithms FastCompSum and CompDot should be preferably used for a computation as with twice the working precision. It can be observed in Fig.4and 5 that, if the condition number is less than about 10 16(K-1) , algorithms SumK and DotK produce results with the maximum possible accuracy in double precision. Then, the accuracy decreases if the condition number increases from about 10 16(K-1) to 10 16K . For condition numbers greater than about 10 16K , results computed by SumK and DotK have no more correct digits. Figures4 and 5are consistent with the properties of algorithms SumK and DotK given in Sect. 7 and 8 for directed rounding: results are computed as in K-fold precision.

Figure 4 :

 4 Figure 4: Accuracy estimated by CADNA using the Sum and the SumK algorithms with 200 randomly generated floating-point numbers.

Figure 5 :

 5 Figure 5: Accuracy estimated by CADNA using the Dot and the DotK algorithms to compute the dot product of arrays of 100 randomly generated elements.

  With rounding to nearest, if nu < 1, the result res of Algorithm 1 satisfies

					(4.3)
	With directed rounding, if nu < 1 2 , then
	|res -s| ≤ γ n-1 (2u)S.	(4.4)
	In Corollary 4.2 Equations 4.3 and 4.4 are rewritten in terms of the condition number on	p i :
	cond	p i =	S |s|	.
	Corollary 4.2.

1

URL address: http://cadna.lip6.fr

Proposition 4.1 ([11]

). Let us suppose Algorithm 1 is applied to floating-point numbers p i ∈ F, 1 ≤ i ≤ n. Let s := p i and S := |p i |.

With rounding to nearest, if nu < 1, then |res -s| ≤ γ n-1 (u)S.

  7: res ← π n + σ n Algorithm 3: Compensated summation of n floating-point numbers p = {p i } using FastTwoSum Proposition 4.3 ([20]). Let us suppose Algorithm 3 (FastCompSum) is applied, with rounding to nearest, to floating-point numbers p i ∈ F, 1 ≤ i ≤ n. Let s := p i and S := |p i |. If nu < 1, then, also in the presence of underflow,

	|res -s| ≤ u|s| + γ 2 n-1 (u)S.
	Corollary 4.4 ([

  Let us suppose Algorithm 3 (FastCompSum) is applied, with directed rounding, to floating-point numbers p i ∈ F, 1 ≤ i ≤ n. For i = 2, ..., n, let e i be the error on the floating-point addition of π i-1 and p i :π i + e i = π i-1 + p i . If nu < 1 2 , then

	Equation 4.5 in Lemma 4.6, established in [8], is recalled for later use.
	Lemma 4.6 ([8]). n	n
			|e i | ≤ γ n-1 (2u)	|p i |.	(4.5)
			i=2	i=1
	A bound for the absolute error on the result of Algorithm 3 (FastCompSum) obtained with
	directed rounding is recalled in Proposition 4.7.
	Proposition 4.7 ([8]). Let us suppose Algorithm 3 (FastCompSum) is applied, with directed round-
	ing, to floating-point numbers p i ∈ F, 1 ≤ i ≤ n. Let s :=	p i and S :=	|p i |. If nu < 1 2 , then,
	also in the presence of underflow,
		|res -s| ≤ 2u|s| + 2(1 + 2u)γ 2 n (2u)S.
	From Proposition 4.7, a bound for the relative error on the result of Algorithm 3 (FastCompSum)
	obtained with directed rounding is deduced in Corollary 4.8.
	Corollary 4.8. With directed rounding, if nu < 1 2 , then, also in the presence of underflow, the
	result res of Algorithm 3 (FastCompSum) satisfies
		|res -s| |s|	≤ 2u + 2(1 + 2u)γ 2 n (2u) cond	p i .
	From Corollary 4.8, because γ n (2u) ≈ 2nu, the relative error bound is essentially (nu) 2 times
	the condition number plus the inevitable rounding 2u due to the working precision.
	The impact of a directed rounding mode on the compensated summation based on the PriestTwoSum
	algorithm is analysed here. The error bounds obtained will be used in Sect. 7. The PriestTwoSum
	algorithm [21, p.14-15] is recalled as Algorithm 4.
	function [c, d] = PriestTwoSum(a, b)
	1: if |b| > |a| then
	2:	exchange a and b
	3:

  Theorem 6.1. Consider a polynomial p of degree n with floating-point coefficients, and a floatingpoint value x. With rounding to nearest, the forward error in the compensated Horner algorithm is

	such that
	|CompHorner(p, x) -p(x)| ≤ u|p(x)| + γ 2 2n (u) p(x).	(6.31)
	It is interesting to interpret the previous theorem in terms of the condition number of the
	evaluation of p at x. Combining the error bound (6.31) with the condition number (6.30) for
	polynomial evaluation gives
	|CompHorner(p, x) -p(x)| |p(x)|	≤ u + γ 2 2n (u) cond(p, x).	(6.32)

  and a bound for the relative error in Corollary 8.2. Proposition 8.1 ([20]). Let floating-point numbers x i , y i ∈ F, 1 ≤ i ≤ n, be given and assume 8nu ≤ 1. Denote by res ∈ F the result computed by Algorithm 13 (DotK) with rounding to nearest. Then |res -x T y| ≤ u + 2γ 2 4n-2 (u) |x T y| + γ K 4n-2 (u)|x T ||y|. Corollary 8.2 ([20]). Assume 8nu ≤ 1. The result res of Algorithm 13 (DotK) obtained with rounding to nearest satisfies res

  | ≤ |x T y| + γ 2n (2u)|x T ||y|.

	From Equations 8.48 and 8.49, we have
	2n
	|r i
	i=1
	(2u) 1 -2n(2u)	|x T ||y|,
	≤ γ 2n (2u)|x T ||y|.	(8.49)

  2u)|x T ||y|. From Proposition 8.3, a bound for the relative error on the result of Algorithm DotK obtained with directed rounding is deduced in Corollary 8.4. Corollary 8.4. Assume 16nu ≤ 1. The result res of Algorithm 13 (DotK) obtained with directed rounding satisfies

	res -x T y x T y	≤ 2u + 2γ 2 4n-2 (2u) +	1 2	γ K 4n-2 (2u) cond(x T y).

  are consistent with the properties of algorithms SumK and DotK given in Sect. 7 and 8 for directed rounding: results are computed as in K-fold precision.

		16
	number of exact significant digits	2 4 6 8 10 12 14					Sum SumK, K=2 SumK, K=3 SumK, K=4 SumK, K=5 SumK, K=6 SumK, K=7
		0
		10 0	10 20	10 40	10 60	10 80	10 100	10 120
					condition number

120 number of exact significant digits condition number Dot DotK, K=2 DotK, K=3 DotK, K=4 DotK, K=5 DotK, K=6 DotK, K=7

Acknowledgement

This work was partly supported by the project FastRelax ANR-14-CE25-0018-01. The authors also wish to thank EDF (Electricité De France) for its financial support.

in terms of execution time varies from 5 to 17 if self-validation is activated. This overhead is higher if any instability is detected mainly because of the heavy cost of cancellation detection.


Lei Chen
email: lei.chen@univ-ubs.fr

Sylvie Gibet
email: sylvie.gibet@irisa.fr

Pierre-Franc ¸ois Marteau

Fabrice Marandola
email: fabrice.marandola@mcgill.ca

Marcelo M Wanderley
email: marcelo.wanderley@mcgill.ca

Quantitative Evaluation of Percussive Gestures by Ranking Trainees versus Teacher

Keywords: percussive gestures, classification, expressive variations, evaluation. ACM Classification H.5.5. Information Interfaces and Presentation (e.g., HCI): Sound and Music Computing: Signal analysis, synthesis, and processing

In this paper we characterize timpani gestures by temporal kinematic features, containing most information responsible for the sound-producing actions. In order to evaluate the feature sets, a classification approach is conducted under three main attack categories (legato, accent and vertical accent) and sub-categories (dynamics, striking position). Two studies are carried out: intra-subject and inter-subjects classification. Results are presented in terms of a quantitative ranking of students, using professional gestures as training set, and their gestures as test set.

INTRODUCTION

While performing a music instrument, musicians establish a rich, well-structured, expressive interaction with the instrument. In order to master such gestural interaction and to control the fine-tuning of the sound-producing gestures, many training years are necessary. In such a learning process, a strong coupling between sound effects and control gestures is established. Instrumental gestures are progressively refined Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org. MOCO'16, July 05-07, 2016, Thessaloniki, GA, Greece c 2016 ACM. ISBN 978-1-4503-4307-7/16/07...$15.00 DOI: http://dx.doi.org/10.1145/2948910.2948934 so that the produced sounds satisfy the desired goals. This process is usually guided by a professional teacher who gives the key elements to improve the regularity and the quality of the performance. In most of instrumental playing, the study of pedagogy around musical training has led to a gestural categorization that follows the gesture-sound relationships.

In previous work a wide range of musical gestures has been studied through the analysis of the gestural signals that are responsible for the sound production, but not so much concern percussive gestures [START_REF] Dahl | On the beat: Human Movement and Timing in the Production and Perception of Music[END_REF], [START_REF] Bouënard | Synthesis of Music Performances: Virtual Character Animation as a Controller of Sound Synthesis[END_REF], [START_REF] Dahl | Striking movements: A Survey of Motion Analysis of Percussionists[END_REF]. Yet these gestures are interesting because they are short striking gestures which can be largely varied along the musical scores, depending on the musical intention. For percussive gestures, different playing modes are usually defined according to the following axes/dimensions: (i) Attack; (ii) Dynamics; (iii) Striking position. Attack refers to the initial stage of sound envelope, which is directly related to sound quality, especially in percussion. In music, expressive terms such as legato, tenuto, staccato, etc. correspond to different types of attacks and indicate various ways of shaping a sound. Dynamics (sometimes called intensity) corresponds to the loudness of the produced sound. The striking position indicates the location where the stroke is played. It is interesting to note that these dimensions may characterize gesture or sound descriptors, or both of them, which is due to the strong gesture-sound coupling that drives percussive gestures. We will also consider this categorization along the above axes and study the effects of distinct expressive variations on the performance.

In this paper we focus on timpani gestures. These gestures are characterized by large kinematic and expressive variations, with all the upper-body articulations of the percussionist strongly involved during gesture executions. Furthermore, the whole kinematics of the timpani gesture embeds not only the sound-producing gesture, but also the preparing and retracting gestures that are supposedly used to ease the performance while reducing the produced effort.

Prior knowledge in musicology and percussion teaching has led to the identification of different sets of gestural attributes in articulation and phrasing for musical performance on timpani. Cook [START_REF] Cook | Teaching percussion[END_REF] proposes the following attributes: (i) manipulating the height of the stroke, which naturally results in grip alterations in order to change the stroke height while keeping dynamic balance; (ii) manipulating the grip, which involves varying finger contact and grip pressure on the stick and results in changes of the stiffness of the system made by the combination of the mallet and the hand of the player; (iii) altering the place where the strokes hit the surface, which involves modifying the vibrations of the playing area; (iv) manipulating the lift of the stroke by changing the arm, wrist and finger movements (configuration and biomechanical state), thus altering the stroke's velocity and time of contact of the head of the mallet with the head of the timpani. F. Marandola proposes to consider: (i) the tightness of the grip depending on the context. This tension has a significant impact on the fullness of the sound, which can be related to the richness of the overall timbre. Hence not enough grip results in less fundamental, and too much grip results in stiffer attack, and greater risk of sound distortion, as well as a reduced range of expressivity; (ii) the depth of the stroke, which is directly linked to the contact time of the stick with the timpani?s skin Inspired by these works, we propose to describe percussive gestures by kinematic temporal trajectories that best characterize these gestural attributes from motion captured data. Figure 1 shows a phase-diagram of the velocity versus position of the right hand mallet 3D-trajectory when a performer plays repeated legato, accent and vertical accent attack beats. We can observe the different patterns achieved when executing different types of repeated attacks. The first purpose of our study is to find out a minimum subset of features that effectively describe timpani gestures, both spatially and temporally. However, to capture a complete human motion, many body segment joints need to be considered, and to describe the angular values or 3D positions of the joints, many attributes are needed accordingly. It is therefore necessary to reduce the dataset by finding spatial as well as temporal reduced representations, while preserving the expressive qualities of the movements responsible for the variations along the percussive performance. The second purpose of our study is to evaluate different performers on the basis of their kinematic similarities with a professional musician. In particular different students' performance are compared to their professor's ones, resulting in a ranking between them according to their capability to reproduce the teachers gestures. A study on percussive gestures is conducted and validated through a classification approach under three main attack categories (legato, accent and vertical accent) and sub-categories (dynamics, striking position), the first one for intra-subject, the second one for inter-subjects classification.

RELATED WORK

Previous research has focused on the analysis of gestures for particular instrumental gestures. Most of the work is directed toward the design of new musical and interactive interfaces. Especially regarding percussion-related systems, an important research direction is the development of devices to track performer gestures for controlling sound synthesis processes [START_REF] Kapur | The Electronic Tabla Controller[END_REF]. Despite the availability of various devices, the most accurate system for tracking percussive gestures remains camera tracking systems [START_REF] Tindale | A Comparison of Sensor Strategies for Capturing Percussive Gestures[END_REF]. These systems offer an effective method for capturing, analyzing and virtually reconstructing the whole body of a performer. Hence, using motion captured (MoCap) data, it becomes possible to identify the gesture profiles and the gesture characteristics that are responsible for the sound production.

One of the issues addressed by the analysis of gestures is to better apprehend the underlying gestural processes involved in the control of sound, and the way it is related to motion data recorded during real performances. A traditional approach consists in identifying from MoCap data the set of features that best characterize the movement. A state of the art of computable descriptors of human motion is described in [START_REF] Larboulette | A Review of Computable Expressive Descriptors of Human Motion[END_REF], presenting both low-level descriptors that compute quantities directly from raw motion data, and higher level descriptors that qualify the meaning, style or expressiveness of 3D motion capture data.

Dahl [START_REF] Dahl | Striking movements: A Survey of Motion Analysis of Percussionists[END_REF] thoroughly studied the striking of drum performance, and tried to establish a correspondence between movement parameters and sound production. The preparatory movement, the rebound and the timing of a stroke were respectively analyzed to study their control over the sound properties. In [START_REF] Dahl | The playing of an accent-preliminary observations from temporal and kinematic analysis of percussionists[END_REF], percussive movement and timing strategies used by professional percussionists were observed. Despite large differences between the preparatory movement of different subjects, results showed a larger preparation height for the accented strokes. Further work [START_REF] Dahl | Playing the Accent: Comparing Striking Velocity and Timing in Ostinato Rhythm Performed by Four Drummers[END_REF] on striking velocity and timing were conducted with several percussion players performing accented strokes at different dynamic levels, tempi and on different striking surfaces attached to a force plate. A consistent individual movement pattern was maintained for the different players under different playing conditions. Another result showed that an increasing height and striking velocity led to an increasing dynamic level, and to a lengthening of the starting interval with the accented strokes.

Timpani gestures have been studied for the purpose of modeling and animating a physical virtual character capable of interacting with a physical sound synthesis system [START_REF] Bouënard | Synthesis of Music Performances: Virtual Character Animation as a Controller of Sound Synthesis[END_REF], [START_REF] Bouënard | Hybrid inverse motion control for virtual characters interacting with sound synthesis -Application to percussion motion[END_REF]. The timpani gesture sequences were cut into elementary beat units, and edited so that a new score could be used to animate the virtual character. Two evaluations were conducted. A qualitative evaluation of synthesized timpani performances, following instructed exercises [START_REF] Bouënard | Virtual Control and Synthesis of Music Performances: Qualitative Evaluation of Synthesized Timpani Exercises[END_REF]. For the quantitative evaluation of the synthesis system, a classification approach was developed to recognize the different attacks under different expressive variations [START_REF] Bouënard | Virtual Gesture Control of Sound Synthesis: Analysis and Classification of Percussion Gestures[END_REF]. As an intuitive hypothesis states that percussionists more specifically control the motion of mallets over time, mallet extremity trajectories were used to conduct the analysis, and more particularly discrete local extrema extracted from position trajectories during beat-to-beat phases [START_REF] Bouënard | Analysis of Percussion Grip for Physically Based Character Animation[END_REF].

A more recent study was conducted [START_REF] Cuykendall | Hearing movement: How taiko can inform automatic recognition of expressive movement qualities[END_REF] on Taiko, a Japanese form of ensemble drumming. Various categories of expressive Taiko performances were captured, recording both drumstrokes (position of the wrist) and sounds. Using machine learning methods, the authors have classified key aspects of Taiko technique, and showed that gesture and sound classification share similar results.

In this paper, we propose to refine our previous work on timpani gestures, and to compare various features of skeletal motion to analyze timpani performances, without any a priori on the kinematic features. Hence we evaluate through a classification approach different sets of joints, of kinematic measures (3D position, velocity, acceleration, etc.), and also how the starting and ending of the beat influences the performance. Furthermore, we propose a methodology to automatically evaluate the different students, using professional gestures as training set, and students gestures as test sets, with a ranking as output.

TIMPANI GESTURES DATABASES AND MOTION REPRE-

SENTATION

Within the three main classes of timpani gestures categorization, i.e. -Attack, Dynamics, Striking position -, we consider the following expressive variations.

• Attack: three playing modes in attack are considered in our study: legato(l), accent(a) and vertical accent(v).

In music performance and notation, legato indicates that musical notes are played smoothly and connected. An accent adds emphasis to a particular note, requiring that it should be played louder than unaccented notes at the same dynamic level. A vertical accent is considered as accented tenuto, which means to hold the note along its full length and play it slightly louder. In terms of gestures, a performer plays legato accent and vertical accent with different positions, velocities and accelerations of the mallet extremity.

• Dynamics: dynamic indications in music are graduated between p or piano, meaning soft, and f or f orte, meaning loud. More subtle degrees of loudness or softness are indicated by: mp, standing for mezzo-piano, meaning moderately soft ; mf , standing for mezzo-f orte, meaning mod-erately loud. The classes of dynamic levels covered in our study are p, mf , and f .

• Striking positions: the basic playing area on a timpani head is at least 4 to 5 cm from the edge of the bowl directly in front of the player, which is nearly 1/3 to the rim. In our study, we request the participants to hit the timpani at the center location (c), at 1/3 to the rim (1/3), or at rim (r) of the timpani. For each class and each variation there is one sequence of 30 strokes. For both databases the performances were executed at the same tempo of 80 bpm, and a pause occurred between the sequences.

Motion Representation

Each recorded sequence is composed of repeated beat units which are manually segmented, each beat being centered on the striking position. This means that one beat contains several frames before the stroke, corresponding to the preparatory movement, and several frames after the stroke, corresponding to the retractive movement. Each beat is a multi-dimensional time series, which is represented as a sequence of k upper-body postures including the stick configuration (with k around 350). Each posture,

X i = {x 1 , x 2 , ..., x m },

Beat Descriptors

Once the beats are segmented, they can be represented by spatio-temporal descriptors. However, as a beat is represented by a matrix of k postures X i , i = 1..k, it is necessary to find out a reduced-dimensional representation that still contains the expressive kinematics and its variation along the beat. In this paper, inspired by knowledge in musicology and percussion teaching, we will assume that the position and velocity of the end position of the mallet in the vertical plane contain sufficient information to characterize timpani performances.

In order to validate our beat descriptors, we use a classification approach. This classification is carried out under various features sets: (i) the position X = (x, y, z) of the mallet extremity; (ii) the position and velocity (dx/dt, dy/dt, dz/dt) of the mallet extremity; (iii) the position, velocity, and acceleration (d 2 x/dt 2 , d 2 y/dt 2 , d 2 z/dt 2 ) of the mallet extremity.

We will also focus on the temporal aspects of the beat, i.e.

we will analyze which part of the beat before and after the striking position mostly influences the classification results. Figure 2 illustrates the z-position trajectories of several superimposed beats achieved for different attack conditions for the referent S 0 (top) and the trainees (S 1 to S 6 ). In this figure we can observe the regularity of the different strokes executed for a given attack, and the changing position of the velocity peaks according to the attack. In the next section we will analyze the intra-subject classification of the main variations of attack (both for the referent S 0 and trainees S 1 to S 6 ). We will also analyze for inter-subject classification the influence of the center position of the beat, and the influence of the window's length taken from either side of the striking position.

Classification

We adopt a simple non-parametric classification method, i.e. As the timpani beats are produced at the same tempo, and are all segmented around the striking moment, they are temporarily aligned and have the same length. Therefore, to evaluate the similarity between the beats, a simple Euclidian distance can be applied on the time series representing them. As experimentally verified, more sophisticated elastic similarity measures such as Dynamic Time Warping do not perform better than the straightforward Euclidean distance.

As results of the classification tasks, we will obtain confusion matrices, which show the success rate of the k-N N classifier for each class.

INTRA-SUBJECT EVALUATION

The aim of the first study is to find an effective subset of motion features that describes the kinematic of timpani gestures and is sufficient to classify the main variations of attack for the referent data S 0 . We first classify the whole set of beats according to the 3 main variations of attack (legato, accent, vertical accent) executed in different conditions of intensity and striking position. The beat descriptors are the z-position of the extremity of the right hand mallet expressed in the shoulder coordinate frame. The k-N N classifier outputs a score rate of 100% for this task.

We then refine this classification task by analyzing the influence of sub-variations related to the intensity (p, mf , f ) or the striking position of the mallet (center, 1/3, rim). The results are showed in Table 1 and 2 with an average classification rate of 98.40% for intensity sub-variations, and 99.83% for striking position sub-variations.

The classification task carried out on subject S 0 applied on the Attack beats (with intensity mf and striking position 1/3) is also validated for the other 6 subjects (S 1 to S 6 ), with the score rates presented in table 3, ranging from 86% to 100%. Hence, the results obtained are very good for intra-subject accuracies. This can probably be explained by the fact that percussionist players have a very long period of learning, resulting in extremely regular skilled gestures.

INTER-SUBJECT EVALUATION

As a result of intra-subject evaluation, we proved that for a large amount of data provided by a performer, the trace of the right hand mallet is able to represent the expressiveness of the gesture. In this section we propose to compare different percussion students in relation to a professional percussionist. The basic idea of inter-subjects evaluation is to use the data of the professional performer S 0 as a training set, and the students data (from S 1 to S 6 ) as a test set, and to classify these data according to attack variations. The evaluation result turns out to be a ranking of how well the students perform comparatively to a professional performer. However, the classification scores fall down to an average of 52.4% when using only the trace of the right hand mallet, as shown in table 4. Therefore, in order to improve the classification, it is necessary to refine the classification experiments.

Following studies on percussion playing [START_REF] Cook | Teaching percussion[END_REF], we consider as descriptors not only the trace of the trajectory of the mallet extremity over time, but also time series expressing the kinematics of the movement (velocity, acceleration, etc.). Hereinafter we will experiment two feature sets, one with the mallet extremity position (x,y,z), and one with both position and velocity of the mallet extremity, i.e. the 6 dimensional time series (x(t),y(t),z(t),vx(t),vy(t),vz(t)).

Furthermore, as observed previously, the beat unit of each subject contains not only sound-producing gesture but also ancillary gesture, which stylizes each performer and tends to introduce a large variability in the performance. Therefore we apply a time window to every stroke, so as to restrict the temporal size of the beat, and thus to focus on the soundproducing gesture while minimizing the effect of the ancillary gesture.

Influence of the Window's Size

In this experiment we consider that the classification is applied on a window centered on the frame when the mallet hits the skin, i.e. frame 175 (middle of the 350-size beats), and we change the size of the window (from 50 frames to 350 frames around the Center 175). The experiment is first conducted with the position of the mallet extremity as feature set, and then with the position and velocity of the mallet extremity.

In the first case (position only), the classification rates are given in table 5 for each subject and window's size, as well as the average for all subjects. We observe that the best results are obtained for a window's size of 150 or 200 for most students (with an average classification rate between 65.8% and 66.6%).

In the second case (position and velocity), the classification scores are given in table 6. We observe that they are sensibly improved (from an average rate of 65.8% for position only to 71.9% for both position and velocity). This supports the hypothesis of Cook, wherein the velocity of the mallet plays an important role in the control of percussive gestures [START_REF] Cook | Teaching percussion[END_REF].

Influence of the Window's Center

To study the importance of the preparatory gesture versus the retractive one, we also tested dissymmetrical windows, i.e.

we analyzed the influence of the center of the window. We tested our classification rates for both feature sets (position and position + velocity), while varying the window's center (from 100 to 250).

The results are respectively given in tables 7 and 8. When selecting position as feature set, the best results are obtained for a window centered around frame 200 (i.e. during the mallet rebound), whereas they are obtained for a window centered around frame 175 for position and velocity (i.e. when the mallet hits the skin). They are also slightly better for the second feature set (72% against 71%).

Besides the improvement of the scores, we also observe more stability within the subjects when selecting the second feature set (position + velocity). Note that we also experimented another feature set comprising position, velocity and acceleration, but the results were about the same as the ones with position and velocity. As acceleration did not improve classification accuracy, we therefore kept this last feature set in the remaining experiments.

EVALUATION THROUGH COMPARISON WITH GROUND TRUTH RANKING

Perceptual evaluations have been conducted by Prof. Marandola to provide a ground truth on the quality of the timpani's performances. This evaluation was carried out on the basis of the videos of the database DB13 showing sequences of beats executed with intensity mf and striking position 1/3. The evaluation focused on both quality of the produced sounds and gestures performed by each student, with a score rating the following sound/gesture attributes (score from 1 to 10).

For the sound was rated 1.1) the fulness of sound; 1.2) the quality of the attack (does the sound correspond to the expected attack?). For the gesture was rated 2.1) the balance between Left and Right Hand; 2.2) the regularity of the dynamics of the strokes; 2.3) the muscular relaxation (absence

Subjects S1 S2 S3 S4 S5 S6
Ground truth rate Gesture Ground truth rate Gesture + Sound Classification rate Position + Velocity of tensions in the neck, shoulders, arms, etc.). For each student separately, and for each category of attack, the different attributes were rated on a scale of 1-10. The results are given in Table 10. The average rates are also given per subject on the basis of the gestural attributes, and the sounds and gestures attributes. We also output a general ranking associated to the average scores.

A comparison between ground truth and automatic classification ranking is then possible. As seen in This experiment shows that focusing the analysis of the trajectory of the mallets on a small window centered on the stroke moment improves the correlation between the classification rates and the perceptual evaluation of the quality of the gestures produced by a subject. One can argue that since this analysis windows corresponds to the most effective part of the sound-producing gesture, it integrates some information that is highly correlated to the sound that is produced and which is used as a perceptual cue by the evaluator.

CONCLUSION

In this paper different sets of raw kinematic features were experimented for describing timpani gestures, and a study about how these features precisely characterize variations of the attack in various contexts was conducted. An automatic classification approach was used for both intra-subject and inter-subjects in varied conditions: different sets of time series (positions, velocities), different window sizes, symmetrical or not around the striking moment.

Our results confirm musical performance hypotheses stressing out the importance of the combination of position and velocity of the mallet trajectories when executing percussion gestures. Note that the 3D trajectories give much better results than only the vertical ones, especially for musicians whose performance is not so regular. In addition, we showed the importance of the window's size for evaluating the quality of the beat attacks: the classification results are more stable and significant around the striking hit (150 frames) than for the whole beat. This result stresses the fact that the most important part of a percussive gesture, during sound production, does not necessarily include ancillary aspects. Ancillary gestures still play an important role, as they ensure that the striking target is reached with the appropriate kinematic, and the wide variation from one performer to another could then be explained by the necessity, for each player, to find the best way to reach the ideal balance between velocity and position, given his/her particular morphology and playing style.

The results are also consistent with the ground truth made from a perceptual evaluation which has led to a ranking of different students on the basis of gestural and sound attributes. This has led us to conclude that raw kinematic trajectories around the striking hit contain most of the information needed for a subjective evaluation. Our approach gives also some insight into a possible quantitative evaluation of students, by classifying their gestures using professional gestures as training set. Therefore, showing that the quality of the students gestures (according to their professor) matches some quantitative score is highly interesting. This could lead to innovative and self-guided ways of learning, that could be adapted to different playing styles pre-recorded from different reference interpreters.

A follow up will be to directly evaluate the individual gestural performance though distance or similarity measures with an expert performer. Concerning the methodology, we intend to use both measures of Precision / Recall to improve the interpretation of the results (threshold effect). We also consider to extend our data set by taking all the beats in various execution situations (varying dynamics and striking positions). Our study could however be extended, taking into account right and left hands, and various performing conditions (variations in intensity, striking position, tempo, etc.). We might also use physical measures about the timing impact with the skin, or compute biomechanical features estimating the relaxation of gesture. Finally, it would be interesting to link this study to a classification of the sound-related signals.

Figure 1 .

 1 Figure 1. Velocity versus position for 3 variations of attack: Legato, Accent, V ertical Accent (data from FM).

Figure 2 .

 2 Figure 2. z-trajectories along time for the subjects S 0 (FM) to S 6 .

  corresponds to the 3D Cartesian positions or the angular values of m = 18 markers located at major joints in the human upper-body plus the position of 3 other markers located at the stick's extremity. Throughout the paper the following notation will be used: each beat B is represented by the triplet (Attack, Dynamics, Striking position), let B = (A, D, S), with A ∈ {l, a, v}, D ∈ {p, mf, f }, and S ∈ {r, 1/3, c}.

Table 1 .

 1 Confusion matrix under 3 attacks and 3 sub-variations of intensity: p, mf , f .

Figure 3 .Figure 4 .

 34 Figure 3. Inter-subject average accuracy according to the subjects, using position and velocity features

  k-Nearest Neighbors, also called k-N N , commonly used in pattern recognition. The data set is divided into a training set and a test set (remaining data). Each sample input from the test set is used to predict the class in which this sample falls, by finding the k closest training examples. The output of the classifier is the class that obtains the majority vote of its neighbors. The nearest neighbors are computed according to a distance between the test beat and the example beats from the training set.

Table 2 .

 2 Confusion matrix under 3 attacks and 3 sub-variations of striking position: center, 1/3, rim.

	Subjects	S1	S2	S3	S4	S5	S6
	Rate	98.11 100 95.35 93.75 86.95 97.87

Table 3 .

 3 Score rates for intra-subject classification for Attack beats (mf , 1/3).

	Subjects	S1	S2	S3	S4	S5	S6
	Rate	41.51 67.35 88.37 35.42 39.13 42.55

Table 4 .

 4 Score rates for inter-subject classification for whole Attack beats.

Table 9 .

 9 Intra-subject ranking: ground truth and automatic classification.

	3	1	2	6	5	4

Table 10 .

 10 Perceptual evaluation: scores (from 1 to 10) using gestural and sound attributes for students S 1 to S 6 .

Table 9

 9

	, we ob-

racy of 37.5. The 3 other subjects S 1 , S 5 and S 6 obtain very close scores for the ground truth, respectively (6.27, 6.67 and 6.87), which is about the same for the classification accuracies (71.70, 65.22 and 65.96). According to these results, it can be argued that the raw kinematic features (position +

 8. Score rates for inter-subject classification using position and velocity features for various window centers, with a window size = 150.