Tom Hubrecht
email: tom.hubrecht@ens.fr

Claude-Pierre Jeannerod
email: claude-pierre.jeannerod@inria.fr

Paul Zimmermann
email: paul.zimmermann@inria.fr

Towards a correctly-rounded and fast power function in binary64 arithmetic

Keywords: IEEE 754, double precision, binary64 format, power function, correct rounding, efficiency

We design algorithms for the correct rounding of the power function x y in the binary64 IEEE 754 format, for all rounding modes, modulo the knowledge of hardest-toround cases. Our implementation of these algorithms largely outperforms previous correctly-rounded implementations and is not far from the efficiency of current mathematical libraries, which are not correctly-rounded. Still, we expect our algorithms can be further improved for speed. The proofs of correctness are fully detailed in the extended version [9] of this paper, with the goal to enable a formal proof of these algorithms. We hope this work will motivate the next IEEE 754 revision committee to require correct rounding for mathematical functions.

I. INTRODUCTION

The IEEE standard for floating-point arithmetic requires correct rounding for basic arithmetic operations since 1985, but in its latest revision (IEEE 754-2019) [START_REF]IEEE. IEEE Standard for Floating-Point Arithmetic (IEEE Std[END_REF], it still does not require correct rounding for the most common mathematical functions such as log, exp, sin, cos, etc.

Among all these commonly used functions, one of the trickiest ones to implement is the power function (x, y) → x y . One difficulty is that it is a bivariate function, thus it has many more possible inputs for a given IEEE format, and in particular its hardest-to-round cases are not known yet for the binary64 format (formerly called 'double precision' and where the precision is 53). Another difficulty of the power function is that it has many exact and midpoint cases (that is, cases where x y is exactly representable in precision 54). Finally, as the relation x y = b y log b x for b > 1 and x > 0 shows, the power function grows very fast in some regions. This requires that some intermediate steps be performed with a significant amount of extra precision which however must be kept small enough in order to achieve high efficiency in practice.

Current mathematical libraries do not provide a correctly rounded power function for the binary64 format. For this format, the maximal known error goes from 0.523 ulp for the GNU libc to 636 ulps for OpenLibm [START_REF] Innocente | Accuracy of mathematical functions in single, double, extended double and quadruple precision[END_REF].

A. State of the Art

Detailed descriptions of how to implement the power function are given in Markstein's book [START_REF] Markstein | IA-64 and Elementary Functions: Speed and Precision[END_REF]Chapter 12] and Beebe's book [START_REF] Beebe | The Mathematical-Function Computation Handbook -Programming Using the MathCW Portable Software Library[END_REF]Chapter 14]. Markstein suggests using the formula x y = 2 y log 2 x , a choice already advocated in [START_REF] Clark | Self-contained exponentiation[END_REF]. We tried this approach first as it looked quite promising for binary arithmetic, but in our context it turned out to be less efficient than using the relation x y = e y log x , essentially since the Taylor approximations of log x and e x have simpler coefficients. Beebe also suggests evaluating x y by means of 2 y log 2 x for binary arithmetic. His algorithms are more complex, in particular since the fused multiply-add (FMA) instruction is not used for producing exact products. Finally, none of these two chapters discusses ways to ensure correct rounding.

The MathLib library [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF] includes a correctly rounded power function for the binary64 format. However, the algorithm used is not fully detailed, and MathLib is no longer maintained. Some MathLib routines (including the power function) were included in the GNU libc up to release 2.27, but the "slow path" was gradually removed, and the latest release of the GNU libc does no longer provide a correctly rounded power [START_REF] Innocente | Accuracy of mathematical functions in single, double, extended double and quadruple precision[END_REF]. Furthermore, MathLib only provides correct rounding for rounding to nearest-even.

Sun Microsystems developed a library called LIBMCR that included in its latest version (release 0.9, April 2004) seven functions with claimed correct rounding (presumably for rounding to nearest-even): exp, log, pow, atan, sin, cos, tan. LIBMCR still compiles with modern compilers but has both efficiency and correctness issues (see §VI).

CRLIBM [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF] also provides a correctly rounded power function. Like for MathLib, this function is only implemented for rounding to nearest-even. Also, this function has remained tagged as "experimental" and has indeed some correctness issues (see §VI). Compared to MathLib, though, CRLIBM has a better treatment of exact cases and midpoint cases, based on works by Lauter and Lefèvre [START_REF] Lauter | Arrondi correct de fonctions mathématiques. Fonctions univariées et bivariées[END_REF], [START_REF] Lauter | An efficient rounding boundary test for pow(x, y) in double precision[END_REF].

More recently, the RLIBM project [START_REF] Lim | High performance correctly rounded math libraries for 32-bit floating point representations[END_REF], the LLVM C Library [START_REF]The LLVM C Library[END_REF], and the CORE-MATH project [START_REF] Sibidanov | The CORE-MATH Project[END_REF] have proposed implementations for various correctly-rounded functions. However, at the time of writing, RLIBM only provides binary32 functions, LLVM libc only provides the hypot and log10 functions in binary64, and CORE-MATH provides about 20 binary64 functions, but not x y .

B. Our Contributions

In this paper we present a new algorithm for the binary64 power function with correct rounding (modulo the knowledge of hardest-to-round cases) together with a very efficient implementation. Compared to previous work (MathLib, LIBMCR, and CRLIBM), this algorithm is not restricted to rounding to nearest-even, but also works for the directed rounding modes from IEEE 754.

Furthermore (and as far as we can tell), this is the first time an algorithm for binary64 powering with correct rounding is fully detailed and presented together with its rounding error analysis. The algorithm, which makes heavy use of the FMA instruction available on most modern processors, can be readily implemented.

Finally, our C implementation of the algorithm yields a twofold to four-fold speedup with respect to MathLib, CRLIBM, and LIBMCR.

C. Outline

The paper is organized as follows. Section II gives a highlevel description of our three-phase algorithm for computing correctly-rounded powers for the IEEE binary64 format, together with some notation to be used later on. Section III then provides a detailed description of all the building blocks and algorithms used in the first phase of the powering algorithm, together with explicit error bounds (each of which having a detailed proof given in [START_REF] Hubrecht | Towards a correctly-rounded and fast power function in binary64 arithmetic. Extended version of this ARITH[END_REF]). Section IV only provides a coarse description of the second and third phases, since those follow essentially the same approach as the first phase (up to higher precisions, larger approximation polynomials, and more sophisticated range reduction). Section V discusses the correctness of our implementation on known worst cases as well as on random inputs. Section VI compares our implementation to the MathLib, LIBMCR and CRLIBM libraries, and to the incorrectly-rounded GNU libc power function. We conclude and discuss further work in Section VII.

II. HIGH-LEVEL ALGORITHM

Before approximating x y using exp(y log x), the case x ≤ 0 gets a special treatment according to the IEEE 754-2019 specification [10, §9]. The same is done for the cases where x or y is ±∞ or NaN. (See also [12, §F.10.4.5].) Thus, we assume from now on that both x and |y| are in the range [α, Ω], where α := 2 -1074 denotes the smallest positive subnormal binary64 number, and Ω := (1 -2 -53 ) • 2 1024 denotes the largest finite binary64 number. Then, clearly, log x is either zero or in a subrange of the range [-Ω, -2 -1022 ]∪ [2 -1022 , Ω] of normal binary64 numbers.

Our powering algorithm has three phases. A first phase approximates x y using double-double arithmetic ( §III) in a very efficient way; if the rounding test of the first phase fails, then a second phase approximates x y using 128-bit arithmetic ( §IV); if the rounding test of the second phase fails again, then a third phase approximates x y using 256-bit arithmetic ( §IV). Since worst cases for the binary64 power function are still unknown, the third phase contains a final rounding test. In the very unlikely event where this last test fails, an error message is printed, which guarantees that no incorrect rounding is produced.

The three phases use the same algorithm:

• first we compute an approximation of log x, by first rewriting x = 2 e • t with √ 2/2 < t < √ 2. Then we use Tang's algorithm [START_REF] Tang | Table-driven implementation of the logarithm function in IEEE floating-point arithmetic[END_REF]: a table value r approximates 1/t and yields a small value z ≈ r • t -1. Finally we get log x ≈ e log 2 -log r + log(1 + z), where an approximation to log r is read from a table and where log(1 + z) is approximated using a small-degree polynomial;

• then we multiply the approximation of log x by y, obtaining a value called r; • finally we approximate exp(r) by writing r = (k + i/2 n ) log 2+z, where k is an integer, i/2 n is a small rational, and z is small. This yields exp(r

) = 2 k •2 i/2 n •exp(z)
, where an approximation to 2 i/2 n is read from a table and where exp(z) is approximated using a small-degree polynomial.

A. Notation and basic properties

We write Z, N, F to denote the set of integers, nonnegative integers, and finite binary64 numbers, respectively. We write ulp(x) to denote the unit in the last place of x ∈ F, which for x ̸ = 0 gives the weight of the last bit of the significand of x. More generally, for any real number x such that |x| ≤ Ω, we define its ulp as ulp(x) := 2 max(-1022,⌊log 2 |x|⌋)-52 , that is, max(α, 2 ⌊log 2 |x|⌋-52 ), with the convention that ulp(0) = α.

We denote by • the current rounding mode, which can be to nearest with ties to even, toward zero, toward +∞, or toward -∞. Our algorithm is division free and makes heavy use of operations of the form t ← •(xy + z), with a single rounding; in practice, such operations are performed efficiently by fused multiply-add (FMA) instructions. Unless stated explicitly, all variables such as x, y, z, t above represent finite binary64 numbers.

Our error bounds will be deduced in particular from the following accuracy properties of •: when |xy + z| ≤ Ω, the absolute error of t is bounded as

|t -(xy + z)| ≤ ulp(xy + z) ≤ ulp(t);
if in addition |xy + z| ≥ 2 -1022 (normal range), then the relative error satisfies |t/(xy + z) -1| ≤ ulp(xy + z)/|xy + z| and is thus always at most 2 -52 . Finally, in some proofs we shall also exploit the fact that if |xy + z| ≤ 2 e for some integer e ≥ -1021, then |t -

(xy + z)| ≤ 0.5 • ulp(2 e ) instead of ≤ ulp(2 e ).

III. FIRST PHASE

The first phase uses double-double arithmetic for efficiency, and delivers an approximation to x y with about 63 correct bits in the general case (see Algorithm 9). The only special cases are those discussed earlier: the cases y ∈ {0, 1/2, 1, 2} are only checked after the first phase, in order not to slow down the critical path.

In the first phase we make heavy use of the following ExactMul and FastTwoSum algorithms. In the absence of underflow and overflow, ExactMul rewrites the exact product of two finite binary64 numbers as a double-double value, that is, as the unevaluated sum of two finite binary64 numbers. The last instruction ℓ ← •(ab -h) of ExactMul is efficiently

Algorithm 1 (ExactMul) Input: a, b ∈ F Output: h, ℓ such that h + ℓ = ab 1: h ← •(ab) 2: ℓ ← •(ab -h)
implemented using an FMA. If ab ∈ αZ (that is, the exact product ab is an integer multiple of α) and if |ab| ≤ Ω, then h + ℓ = ab exactly. FastTwoSum rewrites the exact sum of two finite binary64 numbers as a double-double value whose first component equals the rounded sum and whose second component is (an approximation to) the associated rounding error.

Algorithm 2 (FastTwoSum) Input: a, b ∈ F with a = 0 or |a| ≥ |b| Output: h, ℓ such that h + ℓ approximates a + b 1: h ← •(a + b) 2: t ← •(h -a) 3: ℓ ← •(b -t)
In the absence of underflow and overflow, it is well known that for rounding to nearest FastTwoSum is an error-free transform, that is, h + ℓ = a + b exactly, while for directed roundings we have in general only an approximation: h + ℓ ≈ a + b. It is shown in [22, Theorem 1] that the corresponding approximation error satisfies |h + ℓ -(a + b)| ≤ 2 -105 |h| and that it is zero when the exponent difference between a and b does not exceed 53. (See also [START_REF] Graillat | Tight interval inclusions with compensated algorithms[END_REF] for a weaker bound.) Also, it is shown in [6, Theorem 3] that the middle computation t ← •(h -a) is always exact, whatever the rounding mode.

Besides FastTwoSum, we shall also use the following Fast-Sum algorithm, which allows us to add an extra error term to the low order term produced by FastTwoSum.

Algorithm 3 (FastSum) Input: a, b h , b ℓ ∈ F with a = 0 or |a| ≥ |b h | Output: h, ℓ such that h + ℓ approximates a + b h + b ℓ 1: h, t ← FastTwoSum(a, b h ) 2: ℓ ← •(t + b ℓ )
Lemma 1: In the absence of underflow and overflow, the pair (h, ℓ) computed by Algorithm FastSum satisfies

h + ℓ -(a + b h + b ℓ ) ≤ 2 -105 |h| + ulp(ℓ).
Proof: The error of FastSum is bounded by the sum of the error of the FastTwoSum call, plus the rounding error in t+b ℓ . The first one is bounded by 2 -105 |h| from [START_REF] Zimmermann | Note on FastTwoSum with Directed Roundings[END_REF], and the second one by ulp(ℓ).

P = z -z 2 /2 + 0x1.5555555555558p-2 z 3 -0x1.0000000000003p-2 z 4 + 0x1.999999981f535p-3 z 5 -0x1.55555553d1eb4p-3 z 6 + 0x1.2494526fd4a06p-3 z 7 -0x1.0001f0c80e8cep-3 z 8
Section III-A computes a double-double approximation h+ℓ of log x, which is multiplied by y in Section III-B to obtain a double-double approximation r h + r ℓ of y log x, and a double-double approximation of exp(r h + r ℓ ) is obtained in Section III-C. Finally Section III-D details how these three steps are used in the first phase.

A. Approximation of log x

To approximate log x for x ∈ [α, Ω] a binary64 number, after some argument reduction described below, we have to approximate log(1 + z) for some binary64 number z such that |z| ≤ 33 • 2 -13 . For this task, we use Algorithm p 1 which computes a double-double approximation p h + p ℓ of log(1 + z) -z, using a degree-8 polynomial (Fig. 1).

Algorithm 4 Algorithm p 1

Input: z ∈ F Output: p h + p ℓ approximating log(1 + z) -z 1: w h , w ℓ ← ExactMul(z, z) 2: t ← •(P 8 z + P 7 ) 3: u ← •(P 6 z + P 5 ) 4: v ← •(P 4 z + P 3 ) 5: u ← •(tw h + u) 6: v ← •(uw h + v) 7: u ← •(vw h ) 8: p h ← -0.5 • w h 9: p ℓ ← •(uz -0.5 • w ℓ ) Lemma 2: Given |z| ≤ 33•2 -13
with z an integer multiple of 2 -61 , the double-double approximation p h +p ℓ to log(1+z)-z returned by Algorithm p 1 satisfies The approximation of log x is done using Algorithm 5. First,

|p h + p ℓ -(log(1 + z) -z)| < 2 -75.492 , with |p h | < 2 -16.
x ∈ F ∩ [α, Ω] is written 2 e • t with e ∈ Z and t ∈ F ∩ (1/ √ 2, √ 2
). Then, following [START_REF] Le | Computing floating-point logarithms with fixed-point operations[END_REF], [START_REF] Tang | Table-driven implementation of the logarithm function in IEEE floating-point arithmetic[END_REF], we reduce the range of t even further by computing z = •(rt -1), where r is a precomputed 9-bit approximation to 1/t obtained from i = ⌊2 8 t⌋, and denoted INVERSE i below. (The value of r is explicitly defined in [9, Appendix A-B]; when i ∈ {255, 256}, we shall simply take r = 1 in order to avoid cancellation when approximating log z -log r.) Note that z can be produced efficiently using an FMA.

Lemma 3:

For any t ∈ F ∩ (1/ √ 2, √ 2), let i = ⌊2 8 t⌋ and r = INVERSE i . Then z = •(rt -1) is exact, |z| ≤ 33 • 2 -13 , and z ∈ 2 -61 Z. Proof: See [9, Appendix A-B].
In Lemma 3, it can be checked that the upper bound 33•2 -13 is optimal if we want z = •(rt -1) to be exact, when using a table indexed by ⌊2 8 t⌋.

In Algorithm log 1, the table value LOGINV i is a doubledouble approximation ℓ 1 +ℓ 2 to -log r, such that ℓ 1 ∈ F is an integer multiple of 2 -42 nearest -log r, and ℓ 2 ∈ F is nearest (-log r)-ℓ 1 . Similarly, LOG2H+LOG2L is a double-double approximation to log 2, to nearest and with LOG2H an integer multiple of 2 -42 .

Algorithm 5 log 1 Input: a binary64 value x ∈ [α, Ω] Output: h + ℓ approximating log x 1: write x = t • 2 e with t ∈ (1/ √ 2, √ 2) and e ∈ Z 2: i ← ⌊2 8 t⌋ ▷ i integer 3: r ← INVERSE i , ℓ 1 , ℓ 2 ← LOGINV i 4: z ← •(rt -1) 5: t h ← •(e LOG2H + ℓ 1 ) 6: t ℓ ← •(e LOG2L + ℓ 2 ) 7: h, ℓ ← FastSum(t h , z, t ℓ ) 8: p h , p ℓ ← p 1(z) 9: h, ℓ ← FastSum(h, p h , •(ℓ + p ℓ )) 10: if e = 0 then h, ℓ ← FastTwoSum(h, ℓ) Lemma 4: Given x ∈ [α, Ω], Algorithm log 1 computes (h, ℓ) such that |ℓ| ≤ 2 -23.89 |h|, and |h + ℓ -log x| ≤ ε log • | log x| (1)
with

ε log = 2 -73.527 if x ̸ ∈ (1/ √ 2, √ 2
), and ε log = 2 -67.0544 otherwise. Proof: If x = 1, then it is easy to see that e = 0, t = 1, i = 256 thus r = 1 (remember for i ∈ {255, 256} we use r = 1 to avoid cancellation), z = 0, so that Algorithm p 1 returns p h = p ℓ = 0 and Algorithm log 1 returns h = ℓ = 0. The proof is thus finished in this case.

Let us now assume x ̸ = 1. Firstly, from the proof of Lemma 3, the value z computed at line 4 fulfills the conditions of Lemma 2. The main analysis of Algorithm log 1 is performed in [9, Appendix A-C], where we distinguish three cases:

1) case e ̸ = 0, that is, x < √ 2/2 or √ 2 < x. This case is detailed in [9, Appendix A-D];

2) case e = 0 and i ̸ = {255, 256}, that is, √ 2/2 < x < 255/256 or 257/256 ≤ x < √ 2. This case is detailed in [9, Appendix A-E]; 3) case e = 0 and i ∈ {255, 256}, that is, 255/256 ≤ x < 257/256. This case is detailed in [9, Appendix A-F].

B. Multiplication by y

Once we have computed a double-double approximation h + ℓ to log x, we multiply it by y in order to obtain an approximation to y log x, as in Algorithm 6 below. As pointed out in [START_REF] Markstein | IA-64 and Elementary Functions: Speed and Precision[END_REF]Chapter 12], what is important is to bound the absolute error in the approximation to y log x.

Algorithm 6 Algorithm mul 1

Input: a double-double value h + ℓ, and a double y Output: r h + r ℓ approximating y(h + ℓ)

|r h + r ℓ -y log x| ≤ ε mul with ε mul = 2 -63.799 if x ̸ ∈ (1/ √ 2, √ 2

C. Final exponentiation

Finally we approximate exp(r h + r ℓ ). After some argument reduction described in Algorithm 8, we have to approximate exp(z) for z a double-double value with |z| < 0.000130273. (This bound comes from the paragraph "About the values z h and z ℓ " in [9, Appendix A-I].) We use a degree-4 polynomial (Fig. 2), that is evaluated using Algorithm q 1. Lemma 6: Given (z h , z ℓ ) such that |z h +z ℓ | < 0.000130273 and |z ℓ | ≤ 2 -42.7260 , Algorithm q 1 returns (q h , q ℓ ) such that

Q = 1 + z + z 2 /2 + 0x1.5555555995d37p-3 z 3 + 0x1.55555558489dcp-5 z 4
q h + q ℓ exp(z h + z ℓ ) -1 < 2 -74.169053
and |ℓ| ≤ 2 -42.7096 . Proof: See [9, Appendix A-H].

Algorithm 7 Algorithm q 1 Input: a double-double value z h + z ℓ Output: q h + q ℓ approximating exp(z h + z ℓ )

1: z ← •(z h + z ℓ ) 2: q ← •(Q 4 z h + Q 3 ) 3: q ← •(qz + Q 2 ) 4: h 0 , ℓ 0 ← FastTwoSum(Q 1 , •(q • z)) 5: h 1 , s ← ExactMul(z h , h 0 ) 6: t ← •(z ℓ h 0 + s) 7: ℓ 1 ← •(z h ℓ 0 + t) 8: q h , q ℓ ← FastSum(Q 0 , h 1 ,
INVLN2 ← 0x1.71547652b82fep+12 9: k ← ⌊•(r h • INVLN2)⌉ ▷ nearest integer 10: LN2H ← 0x1.62e42fefa39efp-13 11: LN2L ← 0x1.abc9e3b39803fp-68 12: k h , k ℓ ← ExactMul(k, LN2H) 13: k ℓ ← •(k • LN2L + k ℓ ) 14: z h , z ℓ ← FastSum(•(r h -k h ), r ℓ , -k ℓ ) 15: write k = e • 2 12 + i 2 • 2 6 + i 1 with 0 ≤ i 2 , i 1 < 2 6 16: let h 2 + ℓ 2 approximate 2 i2/64 17: let h 1 + ℓ 1 approximate 2 i1/2 12 18: p h , s ← ExactMul(h 1 , h 2 ) 19: t ← •(ℓ 1 h 2 + s) 20: p ℓ ← •(h 1 ℓ 2 + t) 21: q h , q ℓ ← q 1(z h , z ℓ ) 22: h, s ← ExactMul(p h , q h ) 23: t ← •(p ℓ q h + s) 24: ℓ ← •(p h q ℓ + t) 25: e h , e ℓ ← •(2 e • h), •(2 e • ℓ) Lemma 7: In the case ρ 1 ≤ r h ≤ ρ 2 , if |r ℓ /r h | < 2 -23.8899
and |r ℓ | < 2 -14.4187 , then the value e h + e ℓ returned by Algorithm exp 1 satisfies 16 .

e h + e ℓ exp(r h + r ℓ ) -1 < 2 -74.
Moreover, |e ℓ /e h | ≤ 2 -41. 7 . Proof: Appendix A-J of [START_REF] Hubrecht | Towards a correctly-rounded and fast power function in binary64 arithmetic. Extended version of this ARITH[END_REF] shows that for ρ 1 ≤ r h ≤ ρ 2 , both e h , 2 e • h and

2 e • (h + ℓ) lie in [2 -991 , Ω].
Appendix A-I of [START_REF] Hubrecht | Towards a correctly-rounded and fast power function in binary64 arithmetic. Extended version of this ARITH[END_REF] analyzes Algorithm exp 1 and proves the bound |e ℓ /e h | ≤ 2 -41.7 , while Appendix A-K of [START_REF] Hubrecht | Towards a correctly-rounded and fast power function in binary64 arithmetic. Extended version of this ARITH[END_REF] establishes the relative error bound 2 -74. 16 for e h + e ℓ .

D. Main algorithm

The main algorithm for the first phase is Algorithm 9 (phase 1), which calls in turn Algorithms log 1, mul 1, and exp 1 seen above. The notation RU (2 -63.797 ) means the rounding upwards of the error bound 2 -63.797 ; this value is precomputed or obtained at compile-time, and any value in F larger than 2 -63.797 works. For the rounding test, we cannot use the algorithms from [START_REF] De Dinechin | On Ziv's rounding test[END_REF] which assume rounding to nearesteven, whereas our goal is to design algorithms working for all rounding modes. Instead, we compute a left bound u and a right bound v with the current rounding mode, and if both u and v round to the same value, x y rounds to that value.

Algorithm 9 Algorithm phase 1

Input: two binary64 numbers x, y with x > 0 Output: the correct rounding of x y , or FAIL

1: ℓ h , ℓ ℓ ← log 1(x) 2: r h , r ℓ ← mul 1(ℓ h , ℓ ℓ , y) 3: e h , e ℓ ← exp 1(r h , r ℓ ) 4: if √ 2/2 < x < √ 2 then ε ← RU(2 -57.579 ) 5: else ε ← RU(2 -63.797 ) 6: u ← •(e h + •(e ℓ -εe h )) 7: v ← •(e h + •(e ℓ + εe h )) 8: if u = v then return u 9: else return FAIL Theorem 1: The value returned by Algorithm phase 1 (if not FAIL) is correctly rounded. Proof: If r h < ρ 1 , see [9, Appendix A-L]. If ρ 2 < r h , see [9, Appendix A-M]. If ρ 1 ≤ r h ≤ ρ 2 , see [9, Appendix A-N].
For the first phase, the following four tables are used, which occupy a total of 6416 bytes:

• a table INVERSE used in Algorithm log 1, with r = INVERSE i a 9-bit approximation of 1/t for i • 2 -8 ≤ t < (i + 1) • 2 -8 , 181 ≤ i ≤ 362, such that rt -1 is exactly representable in binary64. When the interval [i•2 -8 , (i+1)•2 -8 ] contains 1, that is, for i ∈ {255,

IV. SECOND AND THIRD PHASES

The second phase uses 128-bit arithmetic (using two 64bit integer words for the significands), and delivers an approximation with about 113 correct bits. It is very similar to the first phase, except for the computation of log x, where we use a two-step (instead of single-step) argument reduction z = r 1 r 2 t -1 so that |z| ≤ 2 -13 ; then we use a degree-9 polynomial for log(1+z), another two-step argument reduction and a degree-7 polynomial for the computation of exp(r). When the rounding test of the second phase fails, exact and midpoint cases are detected using the algorithm from [START_REF] Lauter | An efficient rounding boundary test for pow(x, y) in double precision[END_REF], where it is shown that all exact and midpoint cases belong to the following set:

S = {(x, y) ∈ F 2 | y ∈ N, 2 ≤ y ≤ 35} ∪ {(m, 2 F n) ∈ F 2 | F ∈ Z, -5 ≤ F < 0, n ∈ 2N + 1, 3 ≤ n ≤ 35, m ∈ 2N + 1}.
We have computed hard-to-round cases for all pairs (x, y) from S. Our relative error bound for the second phase is slightly worse than the bound 2 -117 from [14, Algorithm 1], but no remaining worst case at that point defeats that algorithm. (Note that for the new IEEE mode 'to nearest with ties to away', our algorithms should be easy to adapt, since the only difference with 'to nearest-even' is how to round midpoint cases, thus it suffices to adapt the routine handling them.)

The third phase uses 256-bit arithmetic (using four 64-bit integer words for the significands), and delivers an approximation with about 240 correct bits. For the computation of log x, we use the same two-step argument reduction as in the second phase, with a degree-18 polynomial for log(1+z), and another two-step argument reduction with a degree-14 polynomial for the computation of exp(r). Figure 3 summarizes the parameters of the three phases. In particular, the error bound of 2 -63.797 of the first phase when x / ∈ (1/ √ 2, √ 2) corresponds to a probability of about 1/1000 of calling the second phase, thus to an average overhead of less than one cycle with respect to the first phase and corresponding rounding test.

Lauter estimates to 2 112 the number of pairs (x, y) such that x y lies in the binary64 range [START_REF] Lauter | Arrondi correct de fonctions mathématiques. Fonctions univariées et bivariées[END_REF]. Thus a full search for hardest-to-round cases is not possible, even with clever algorithms like SLZ [START_REF] Stehl É | Searching worst cases of a one-variable function using lattice reduction[END_REF]. As a consequence, it is possible, albeit extremely unlikely, that the rounding test of the third phase fails. In such a case, an error message is printed, with the corresponding inputs x and y, which will enrich the knowledge of hard-to-round cases. This guarantees that whenever the function does not print this error message, the returned value is correctly rounded. If worst-case information was available, one could avoid the rounding test of the third phase, by adding some exceptional cases if needed. This would ensure the function never yields an error message and always returns a correctly-rounded result, but this would have little impact on its efficiency.

V. CORRECTNESS

We have tested the correctness of our implementation, which is publicly available at https://gitlab.inria.fr/zimmerma/ core-math-power-b64, for all rounding modes, on a set of 917231 input pairs (x, y). This set includes: (a) worst cases for exponents y that might yield exact or midpoint cases, including cases where x y lies in the subnormal range; (b) special values specified by IEEE 754-2019, for example (x, y) = (1, -0) or (±∞, -0); (c) pairs with x < 0 and y an integer; (d) pairs with x near 1 and x y near underflow/overflow; (e) exact and midpoint cases; (f) inputs exhibiting bugs in other libraries; (g) non-regression tests; (h) additional worst cases found in the literature, in particular from Section 14.14 in [START_REF] Beebe | The Mathematical-Function Computation Handbook -Programming Using the MathCW Portable Software Library[END_REF], and worst cases computed by Vincent Lefèvre for x n and x 1/n , n integer. We also tested our implementation on a set of 10 8 random inputs pairs, again for all rounding modes. For all these tests, we used as reference the value given by GNU MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF], with exponent range matching that of binary64, and emulating subnormals using mpfr_subnormalize.

VI. EFFICIENCY

In this section, we measure the efficiency of our implementation in the C language of the algorithms described in this paper. We measure both the reciprocal throughput and the latency, using the CORE-MATH perf.sh tool [START_REF] Sibidanov | The CORE-MATH Project[END_REF]. For the power function, CORE-MATH generates random x and y uniformly in [0, 20], so that x y lies in [0, 10 26 ] approximately. We compare our implementation to the GNU libc (which is not correctly rounded), to MathLib, CRLIBM, and CRLIBM, for rounding to nearest-even.

For MathLib, we used the non-official copy from https: //github.com/dreal-deps/mathlib, which we compiled with the default compile flags (using -O3 -march=native yields plenty of incorrect roundings). Our tests confirm that the MathLib upow routine does not provide correct rounding for directed rounding modes.

For LIBMCR, we used the non-official copy from https: //github.com/simonbyrne/libmcr/, which matches our copy of release 0.9 from 2004. For rounding to nearest-even, on our set of about one million tests, the pow function from LIBMCR does not terminate for 15 inputs pairs, for example x = 0x1.470574d68e0afp+1 and y = 0x1.02e0706205c0ep+1, and we get about 96% of failures on the remaining tests, for example for x = 0x1.f80b060553772p-1 and y = 0x1.99cp+13, LIBMCR yields 0x1.00001p+0. The power function from CRLIBM is explicitly said as experimental, and only works for rounding to nearesteven. We also noticed some issues, for example for x = 0x1.524ebae943097p+1 and y = 0x1.ep-2, it returns -5 instead of 0x1.93bd0cd47eb5fp+0. Looking at the code, it appears that the strange return value of -5 corresponds to a failure in the last rounding test, which is not treated (here, x y has 68 identical bits after the round bit).

GNU libc

Figure 4 shows that our implementation is only about 50% slower than the GNU libc, which is not correctly rounded (about 40% for the latency). (On a i7-1260P, we get a reciprocal throughput of 35 cycles and a latency of 78 cycles, against 23 and 59 respectively for the GNU libc.) While Figure 4 only considers rounding to nearest-even, timings of our implementation are very similar for other rounding modes. For the reciprocal throughput, our implementation outperforms that of MathLib by 86%, that of CRLIBM by a factor larger than 3, and that of LIBMCR by a factor of almost 4.

VII. CONCLUSION

We have proposed algorithms for computing the power function x y with correct rounding, for the binary64 floatingpoint format and the four main rounding modes (to nearesteven as well as the three directed roundings of IEEE 754), and up to the knowledge of hardest-to-round cases. Our approach makes heavy use of the FMA instruction (either in hardware or emulated in software), and leads to a very efficient implementation, which is 2x to 4x faster than the power functions from MathLib, CRLIBM, and LIBMCR, claiming correct rounding.

When designing these algorithms, our aim has been to ensure both correctness and high efficiency. Efficiency follows from a three-phase approach whose first phase handles most binary64 inputs (x, y) very fast thanks to carefully designed FMA-based double-double computations. Correctness is obtained by returning either the correctly-rounded value of x y (for which detailed error analysis has been done for the first phase), or switch to the next (more accurate) phase. Our implementation handles all the special input cases (such as x ≤ 0, x or y NaN or infinite, etc.) according to the current IEEE 754 specification [10, §9.2.1].

As can be seen from the material in [START_REF] Hubrecht | Towards a correctly-rounded and fast power function in binary64 arithmetic. Extended version of this ARITH[END_REF], we have systematically based our algorithms and implementations on very detailed proofs, which allows us to carefully control the accuracy of each of the routines involved as well as to predict and handle the (im)possibility of underflow and overflow. Furthermore, in the unlikely case where we cannot conclude that the correct value is returned for x y , returning instead an error message along with the corresponding input (x, y) makes it possible to usefully augment the list of known hardest-toround cases for that function.

Although our current implementation is already very fast, we do not claim that it cannot be improved further for speed. We hope that our results will contribute to turn the current IEEE recommendation of correct rounding for the mathematical functions listed in [10, Table 9.1] into a requirement. Our approach also easily extends to other rounding modes (such as round to nearest with ties to away) and to other formats (such as double extended and binary128, or decimal formats). Finally, the high level of detail of our proofs will hopefully make it easier to develop formal proofs of our algorithms, and we are planning to investigate this in the near future as well.

Fig. 1 .

 1 Fig.1. Degree-8 polynomial P (z) generated by Sollya[START_REF] Chevillard | Sollya: an environment for the development of numerical codes[END_REF] for approximating log(1 + z) when |z| ≤ 33 • 2 -13 , with absolute error at most 2 -81.63 and relative error at most 2 -72.423 .

9

 9 and |p ℓ | < 2 -25.446 . If z ̸ = 0, and assuming further |z| < 32 • 2 -13 , the relative error satisfies z + p h + p ℓ log(1 + z) -1 < 2 -67.441 . Proof: See [9, Appendix A-A], and refinement via |δ 0 |+|δ 6 | ≤ 3.505u for the relative error in the case |z| < 32 • 2 -13 .

  ), and ε mul = 2 -57.580 otherwise. Proof: See [9, Appendix A-G].

Fig. 2 .

 2 Fig. 2. Degree-4 polynomial Q(z) generated by Sollya for approximating exp(z) when |z| ≤ 0.000130273, with absolute error at most 2 -74.346 .

  If x, h, ℓ are as in Lemma 4 and if the exact product yh satisfies 2 -969 ≤ |yh| ≤ 709.7827, then Algorithm 6 computes (r h , r ℓ ) such that |r h | ∈ [2 -970 , 709.79], |r ℓ | ≤ 2 -14.4187 , |r ℓ /r h | ≤ 2 -23.8899 , |r h + r ℓ | ≤ 709.79, and

1: r h , s ← ExactMul(y, h) 2: r ℓ ← •(yℓ + s) Lemma 5:

  ℓ 1 ) Algorithm 8 Algorithm exp 1 Input: a double-double value r h + r ℓ Output: e h + e ℓ approximating exp(r h + r ℓ ) 1: ρ 0 = -0x1.74910ee4e8a27p+9 ≈ -745.133 2: ρ 1 = -0x1.577453f1799a6p+9 ≈ -686.909 3: ρ 2 = 0x1.62e42e709a95bp+9 ≈ 709.78267 4: ρ 3 = 0x1.62e4316ea5df9p+9 ≈ 709.78276 5: if ρ 3 < r h then return e h = e ℓ = Ω 6: if r h < ρ 0 then return e h = α, e ℓ = -α 7: if r h < ρ 1 or ρ 2 < r h then return e h = e ℓ = NaN 8:

  256}, we use r = 1 to avoid cancellation in the computation of log z -log r. This table has 182 binary64 entries, thus occupies 1456 bytes; • a table LOGINV used in Algorithm log 1 such that for 181 ≤ i ≤ 362, LOGINV i is a double-double approximation to nearest of -log r, where r is defined above. This table has 182 entries with two binary64 values, thus occupies 2912 bytes; • a table T 2 such that for 0 ≤ i < 64, T 2 [i] is a double-double approximation of 2 i/2 6 , used in line 16 of Algorithm exp 1. It has 64 entries with two binary64 values, thus occupies 1024 bytes; • a table T 1 such that for 0 ≤ i < 64, T 1 [i] is a double-double approximation of 2 i/2 12 , used in line 17 of Algorithm exp 1. It has 64 entries with two binary64 values, thus occupies 1024 bytes.

ACKNOWLEDGEMENTS

The authors are grateful to the three anonymous referees for their useful feedback, to Laurence Rideau and Laurent Théry who found some small issues in the initial proof of Lemma 2 while trying to convert it into a formal proof, and to Vincenzo Innocente for his feedback at early stages of this research. The first author is partially funded by CERN, and the search for worst cases was performed using computer resources from CERN and Grid 5000.


Yujiro Okuya
email: yujiro.okuya@limsi.fr

Nicolas Ladeveze
email: nicolas.ladeveze@limsi.fr

Cédric Fleury
email: cedric.fleury@lri.fr

Patrick Bourdot
email: patrick.bourdot@limsi.fr

VR-CAD Framework for Parametric Data Modification with a 3D Shape-based Interaction

Keywords: Virtual Reality, Computer aided design, Parametric CAD, Haptics

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Most of Virtual Reality (VR) applications for Computer Aided Design (CAD) focused on immersive project reviewing or early prototyping in order to reduce cost and time of product design cycle. In terms of design activity, current VR technique enables industrial designers to draw preliminary drafts or edit simple primitives with their hand gestures in 3D space. Such applications enhance designer's inspiration and improve their spatial thinking for aesthetic design. However, hand gesture based interaction has not yet been applied for parametric design tools due to the complexity of CAD data. Industrial models are necessarily modeled with parametric CAD systems (e.g. CATIA, SolidWorks, Inventor, etc.) based on preliminary drafts by CAD engineers, in order to support detailed design requirements and numerical simulations for product's manufacturability. The different focus of industrial designers and CAD engineers impose them to use distinct tools for design activities. This disparity interferes with their communication in early design stage. In this poster, we present a new VR-CAD framework, allowing user to modify parametric CAD data with 3D interaction in an immersive environment. With this framework, users can implicitly modify parameter values of CAD data with co-localized 3D shape-based interaction. This poster describes the system architecture and the interaction technique based on it.

Related Work

Most VR-CAD applications place value on coherent dimension between visualization and interaction space. According to Clark, "To expect a designer of 3-D surfaces to work with 2-D input and viewing devices unnecessarily removes a valuable degree of freedom" [START_REF] Clark | Designing surfaces in 3-d[END_REF]. 3D interaction technique has been mainly carried out for drawing tools and surface modeling, or a few of them are applied for manipulation of 3D primitives. Pushing/pulling interaction for object manipulation was introduced in [START_REF] Dani | Creation of concept shape designs via a virtual reality interface[END_REF], and MockUp Builder [START_REF] De Araújo | Mockup builder: 3D modeling on and above the surface[END_REF] applied this interaction for object creation/edition on multi-touch table top display with co-localized manipulation. However, it is difficult to apply these interaction techniques to complex parametric CAD modeling. A parametric CAD model is defined by a set of operations (Extrusion, Sweep, Boolean, and so on) applied from primitives and 2D sketches, based on a number of parameters and geometrical or topological constraints. In order to update the CAD model, users need to select a specific constraint in Constructive History Graph (CHG) and modify its parameter value. cReaRV [START_REF] Martin | A VR-CAD data model for immersive design: The cRea-VR proof of concept[END_REF] adapted this selection method to 3D space by allowing users to implicitly select a parameter value from selection of a relevant surface with a flystick. Nevertheless, parametric modification still remains in one dimension: users increase or decrease a parameter value by a scroll motion of their hands in 3D space. This one dimensional parametric modification can hardly be transformed into 3D interaction, especially with the complex model containing several constraints.

VR-CAD System Architecture

To address this issue, we designed a new VR-CAD system based on a distributed architecture (Figure 1). It consists of three components: a VR-CAD server, a Visuo-Haptic server (VHServer) and a VR platform. The VR-

14th annual EuroVR conference -EuroVR 2017

Poster track 2 CAD server interacts with a CAD engine (CATIA V5) to load and update any parameter modifications into the CAD data structure. The VHServer handles user interaction during modification (see section 4), and supports force feedback device with a high communication rate. The VR platform manages the graphic rendering of virtual environment (VE) with Unity. Heterogeneous platforms, such as CAVE, Wall-sized display and HMD, are supported. Figure 2 is an example of this VR-CAD application in a CAVE system with a force feedback device.

3D Shape-based Interaction to Modify CAD-data Parameters

Main difficulty of a 3D shape-based interaction with a parametric CAD model is the "unpredictable behavior" of shape evolution due to the internal constraints. Industrial CAD model often contains geometrical or topological constraints between two or more entities to define the design rules or limit its physical position and motion. In order to anticipate the direction of the shape deformation, we compute several meshes by a set of discrete parameter values with a slight offset from the initial value. In such a way, users can choose one of the proposed shape with a 3D hand motion. Only the closest shape from user's hand is rendered in the VE, so that the shape appears to be following user's hand as a pushing/pulling interaction on surface regardless of hidden internal topology complexity. When a user select a parameter by picking a specific part of a CAD object in the VE, the information about the selected part and a set of new parameter values are sent to the VR-CAD server. Then, the VR-CAD server tessellates different B-Rep data computed from the new parameter set with the CAD engine and save generated shapes (.obj files) into a shared folder. The VH server imports mesh data and computes a distance between meshes and user's hand position. It then transmits the closest mesh ID to a VR platform to switch the visualized mesh.

Conclusion and Future Work

This paper presents a new VR-CAD framework for CAD-data modification in an immersive environment. The contributions of this framework are to support (i) native CAD-data modification in VR, (ii) 3D shape-based interaction for parametric CAD-data modification, and (iii) CAD design activity over heterogeneous VR systems. We are currently investigating users' assessment of the proposed interaction technique and collaborative design scenario. We expect that the framework improves the communication between industrial designers and CAD engineers in early design stage by providing a common design space.

Figure 1 .

 1 Figure 1. System architecture of the distributed VR-CAD system.Figure2. Example of the application in a CAVE system with a force feedback device.

Figure 2 .

 2 Figure 1. System architecture of the distributed VR-CAD system.Figure2. Example of the application in a CAVE system with a force feedback device.


Sebastien Roth

Jean-Sebastien Raul

Remy Willinger

Finite element modelling of paediatric head impact: Global validation against experimental data

Biomechanics of the human head has been widely studied for several decades. At a mechanical level, the use of engineering and Finite Element (FE) methods has allowed injury mechanisms to be investigated using biofidelic FE models. These models are generally validated using experimental data then used to simulate real-world head trauma in order to derive numerical tolerance limits, leading to efficient injury predicting tools. Due to ethical issues, experimental tests on the paediatric population remain prohibitive so direct validations of numerical models cannot be performed. However injury biomechanics on paediatric population is emerging with experimental tests on the paediatric cadavers or tests on biological tissue and the development of finite element models. The present paper proposes a new finite element model of a newborn head, simulating its main features, with material properties from the literature. Global validation of the model against experimental data in terms of skull deflection is performed and the model is used to simulate paediatric skull fracture coming from real-world head trauma.

Introduction

Paediatric FE models

According to epidemiological studies [START_REF] Viano | Serious brain injury from traffic related causes: priorities for primary intervention[END_REF], Traumatic Brain Injury (TBI) is a major cause of death and disability among the young paediatric population. Most frequent causes of TBI are vehicle traffic accidents, falls and mishandling (e.g. shaken baby syndrome) [START_REF] Langlois | Traumatic brain injury in the United States : emergency department visits, hospitalizations, and death[END_REF][START_REF] Cdc | Childhood Injuries in the United States[END_REF]. Determination of reliable injury criteria and tolerance limits is essential for the understanding of injury occurrence and for the development of protective devices. Finite element modelling is a widely used tool for investigation of the dynamic response of the head under impact. Compared to adult models, there are relatively few realistic three dimensional finite element models of a child's head.

A finite element model of a 6-month-old child has been developed by DeSantis Klinich et al. [START_REF] Desantis Klinich | Estimating infant head injury criteria and impact response using crash reconstruction and finite element modeling[END_REF] and intracerebral response was analyzed through the simulation of real-world accidents.

The specific injury scenarios in the reconstruction involved infants in rearfacing child restraint systems (CRS) who sustained skull fractures and brain injuries from deploying airbags. These numerical replications allowed developing injury criteria for minor skull fracture. Roth et al. [START_REF] Roth | Finite element analysis of impact and shaking inflicted to a child[END_REF] performed a comparative study between impact and shaking situations with a 6-month-old child head model. This study showed similar results in terms of bridging veins strain for both scenarios. Limitation of scaling methods in child head modelling was performed, comparing a scaled adult model to a biofidelic 3-year-old child model [START_REF] Roth | Limitation of scaling methods in child head finite element modelling[END_REF]. Differences in terms of intracranial parameters led to the conclusion that accurate geometry of child head is needed for a child head numerical investigation. Tolerance limits for neurological lesions were developed with this 3-year-old child head model [START_REF] Roth | Child head injury criteria investigation through numerical simulation of real world trauma[END_REF], by simulating real-world head trauma. Shearing stress appeared to be clearly the most predictive parameters, allowing clear distinction between injured and non-injured cases Coats et al [START_REF] Coats | Parametric study of head impact in the infant[END_REF] developed a 1.5-month-old infant FE model of the head to reproduce Weber's experiments on bone fracture [START_REF] Weber | Experimental studies of skull fracture in infant[END_REF]. Coats et al. compared results of the numerical simulations to experimental tests for different impacted surfaces, in terms of fracture line. Replication of a number of falls allow finding tolerance limit for skull fracture for the newborn infant head. Finally, only DeSantis Klinich and Coats investigated skull fracture with finite element models, by reconstructing several real-world head traumas. However, even if these models give encouraging results in the field of child head biomechanics, none of them have been validated against experimental data.

Experimental tests on paediatric population

Typically, implementation of mechanical properties of head components is based on experimental data. While adult brain and bone material have been extensively investigated, few data are available for the paediatric population: tensile and flexion tests have been performed [START_REF] Margulies | Infant skull and suture properties: measurements and implications for mechanisms of paediatric brain injury[END_REF] in bone and sutures samples extracted from piglets, allowing elaboration of constitutive relationships. Infant bone properties have also been investigated [START_REF] Mc Pherson | The elastic modulus of fetal cranial bone: a first step toward an understanding of the biomechanics of fetal head molding[END_REF][START_REF] Jans | Bending properties of cranial bone segments of new-born children[END_REF][START_REF] Kriewall | Bending properties and ash content of fetal cranial bone[END_REF][START_REF] Kriewall | Structural, mechanical, and material properties of fetal cranial bone[END_REF][START_REF] Runge | Material properties of human infant skull and suture: experiments and numerical analysis[END_REF], showing its orthotropic properties, based on 3-point flexion experiments. The above-mentioned experiments have been performed on quasi-static conditions. High strain rate experiments have also been conducted [START_REF] Coats | Characterization of paediatric porcine skull properties during impact[END_REF][START_REF] Coats | Material properties of human infant skull and suture at high rates[END_REF] and influence of strain rate on elasticity modulus was studied. Experimental tests on brain samples have also been performed [START_REF] Thibault | Age dependent material properties of the porcine cerebrum: effect on paediatric inertial head injury criteria[END_REF]. In parallel to these experimental characterizations on biological tissues, tests on a whole paediatric head have been performed by [START_REF] Prange | Mechanical Properties and Anthropometry of the Human Infant Head[END_REF] in order to investigate the global response of the head submitted to compression tests [START_REF] Prange | Mechanical Properties and Anthropometry of the Human Infant Head[END_REF]. For different impact tests (frontal, occipital, vertex, parietal), impact force and acceleration were computed allowing to give first mechanical results on paediatric cadaver.

Based on the above experimental data, and motivated by the need of a detailed validated newborn head model, the present study proposes the development of a new FE model of newborn child head, based on CT scan slices. In order to ensure its biofidelic behaviour, this model is validated against experimental data and is used in the framework of real-world head accident replication. To the author knowledge, this model is the only one combining biofidelic geometry, global validation against experimental tests, and a use in a context of accident replication, which can contribute to the understanding of injury mechanism in children population.

Method

Geometry

The proposed finite element model was based on the geometrical 3D reconstruction of slices obtained by CT scanner. The specifications of the scanner were not provided. Scanner examination was performed on a newborn subject (17-day-old infant) for medical purposes with no relation with the present study. No abnormalities were found by the medical staff. Resolution of 2D slices was millimetric. For the reconstruction, 230 slices in transversal plane were used. The scanner sections underwent grey level processing in order to distinguish the skull from the soft tissues. A 3D Triangular mesh were then generated and imported into Stereolithographic format and imported in Hypermesh v9.0 software (ALTAIR©) for a step of regular and homogeneous meshing, as needed by explicit codes like Radioss as illustrated in Figure 1.

Anthropometry

Global measurements of the head model were checked based on anthropometric studies concerning the evolution of the head during growth. The main dimensions (length, width, and circumference) of the model were compared to anatomical studies reported by Dekaban et al. [START_REF] Dekaban | Table of cranial and orbital measurements, cranial volume and derived indices in males and females from 7 days to 20 years of age[END_REF] as illustrated in Figure 2. On the other hand, the model's cerebral volume has been compared to anthropometric data (Figure 2), and illustrates that the modelled volume (418 cm 3 ) is realistic. In terms of dimension and mass, the model corresponds to an average population. There are very few discrepancies between the model and the curves, which do not reach more than 3.3 % difference as illustrated in Table 1. The head model can therefore be considered as a 50 percentile newborn head.

Meshing

This model includes the main anatomical features of a newborn head, as described in Figure 1: the scalp, the Cerebro-Spinal Fluid (CSF) and the brain were modelled with brick elements. Coronal, sagittal, lambdoid and temporo-parietal sutures were modelled with shell elements as well as anterior, posterior and sphenoidal fontanels. Skull and membranes were modelled with shell elements. The skull thickness of the model has been evaluated from CT scan slices. In area of the skull which are important in head impact investigations (parietal, frontal, occipital, temporal bones), the thickness is almost constant. Finally, the thickness of the skull is 1.4mm, and is assumed homogeneous, based on measurements given by the CT scan, as it was performed in previous study of head impact [START_REF] Roth | Child head injury criteria investigation through numerical simulation of real world trauma[END_REF].

The whole model is comprised of 30000 solid and 3700 shell elements as shown in Figure 1. The mass of the global head is 0.72 kg.

Material properties

As explained in the introduction, only a few studies report mechanical properties of child head components. Membrane and scalp properties were derived from the literature [START_REF] Zhou | Head injury assessment of a real world crash by finite element modelling[END_REF] concerning adult subjects with a linear elastic constitutive law. Cerebro-Spinal Fluid modelling was considered similar to adult subarachnoid space. Thus mechanical properties reported by Willinger et al were considered in this model [START_REF] Willinger | Modal and temporal analysis of head mathematical models[END_REF]. Coats et al. [START_REF] Coats | Characterization of paediatric porcine skull properties during impact[END_REF][START_REF] Coats | Material properties of human infant skull and suture at high rates[END_REF] have investigated material properties of newborn skull and sutures. The constitutive law of sutures and fontanels were considered as linear elastic based on tension tests. The constitutive law of the skull was elastic-plastic with rupture, based on three-point bending tests. Material properties of the face was considered the same as the skull.

There is a considerable paucity of data concerning the evolution of brain tissue mechanical properties with age. Few studies have investigated properties of paediatric population: Thibault et al. [START_REF] Thibault | Age dependent material properties of the porcine cerebrum: effect on paediatric inertial head injury criteria[END_REF] compared shear viscoelastic properties of 2 to 3-day-old (~equivalent to a 1-month-old human) and 1-year-old (equivalent 4 year-old human) porcine subjects. Non linear aspects of brain tissue have also been investigated in the literature [START_REF] Prange | Regional, Directional, and Age-dependent Properties of the Brain Undergoing Large Deformation[END_REF].

When comparing all these material properties, it appears that some discrepancies exist concerning brain tissues stiffness. In the present study, only linear aspects of brain constitutive law are taken into account. However, in the light of the discrepancies in terms of brain laws, a parametric study is performed on viscoelastic parameters of brain law. Baseline value of material properties implemented in the finite element model are summarised in Table 2.

Mesh convergence

Three different mesh were tested. 5 mm, 2.5 mm (the present model), and a fine mesh ( 1.2 mm). The model meshed with 5mm elements had a quite high time pitch (0.5e-3 ms), but element size may not respect the original geometry of complex head structures. Mesh size was too large. Fine mesh respected the geometry of the brain and the head, but had a low time pitch (8e-5 ms). The present model with element of 2.5mm size, combined both respect of the geometry and acceptable time pitch ( 0.2e-3 ms) for an acceptable computational time. For a 50 ms simulation, run lasted almost 15 hours. Differences in terms of mechanical outputs did not bring significant difference between 1.2 mm and 2.5 mm size elements. Results in terms of stress distribution were similar and only 7% difference were noticed in terms of maximum stress. Simulations performed with the FEM using 2.5 mm size elements showed that the mesh was sufficiently refined to perform a run in a stable manner over the range of loading conditions, without large element distortion, and gave an acceptable computational time. As a result, the mesh size was assumed to be acceptable. Finally, this density was a good compromise between the respect of the geometry and the computational time.

Results

Validation under experimental impact.

Finite element models of the adult head have been typically validated against experimental data carried out in terms of intracerebral pressure, brain displacement and skull fracture [START_REF] Nahum | Intracranial pressure dynamics during head impact[END_REF][START_REF] Hardy | Investigation of head injury mechanisms using neutral density technology and high-speed biplanar X-ray[END_REF][START_REF] Yoganandan | Biomechanics of skull fracture[END_REF]. For paediatric populations, there are very few experimental tests allowing validations of numerical models. To the authors' knowledge, only [START_REF] Prange | Mechanical Properties and Anthropometry of the Human Infant Head[END_REF] [START_REF] Prange | Mechanical Properties and Anthropometry of the Human Infant Head[END_REF] study provides validation data for the global behaviour of the head. In the reported study, the whole newborn paediatric head was positioned between two flat plates and was submitted to a compression in the antero-posterior direction at different velocities along the Frankfort plane. Force and skull deflection were recorded during these tests.

As in the experimental test, the proposed finite element model of the newborn head was subjected to an anterior-posterior compressive force. For all the finite element analysis performed in this study, energy balance was checked and Hourglass energy limit was set to 10% of total energy. Force versus deflection was computed and superimposed to the experimental curve. Figure 3 illustrates the simulation of the compression of the head at 50 mm/s. Numerical and experimental forces versus deflection curves are also reported for the 50 mm/s velocity. Good correlations were found between experimental and numerical results in terms force displacement curves. Prange et al. also performed experimental drop tests with newborn heads, at different velocity on to five different impact locations (Occipital, vertex, forehead, right and left parietal). Impact force versus time was recorded. Dividing impact force by the head mass, acceleration curves were reported in this study. The conclusion of these experimental tests is that average peaks acceleration are significantly different for different drop height but do not vary between impact locations. Then, the 30 cm experimental height drop test was reproduced with the new FE model and led to the computed acceleration versus time curve reported in Figure 4 and superimposed to the experimental test results. Numerical accelerations were also obtained by dividing the impact force by the head mass. Based on these experimental curves coupled with the provided standard deviations, and in order to provide a quantitative model validation, statistical analyses have been performed on the response of the FEM for each impact configuration.

For the forehead impact, a good correlation in terms of acceleration was observed. Indeed, 100% of the response was included in the experimental corridor. For the occipital drop test, as well as for parietal drop test, the error was of the order of 51 %, and 45% respectively, showing a less accurate response. For the vertex drop test, a good correlation was observed between the FEM and the corridor with 85 % of the response included in the experimental corridor. Analysis of simulations for each impact location led to the same conclusion as in Prange's study concerning peak acceleration i.e. does not vary significantly between impact locations.

Parametric study

In the field of biomechanics, large discrepancies exist in the literature in terms brain material properties [START_REF] Nicolle | Shear properties of brain tissue over a frequency range relevant for automotive impact situations: new experimental results[END_REF]. As for adult brain tissues, very few experiments exist for paediatric brain tissues. This point should be investigated in mathematical models. Indeed, it is necessary to know the influence of material law parameters on the model response, as these parameters are not well controlled. Especially, for the present model where skull behaviour is investigated, it is important to study the influence of the viscoelastic parameters of the brain tissue on the skull deformation. Therefore a parametric study was performed on brain characteristics, for a 30 cm drop height forehead impact, as described in Table 2. Parameters of the viscoelastic law G 0 , G ∞ and bulk modulus have been changed in order to investigate their influence on the model's response.

In Table 3, bold values indicate a greater than 15% difference from baseline value A large variation of viscoelastic parameters led to minor variation of mechanical parameters such as Von Mises stress in the skull, skull deflection or peak acceleration of the head. Only higher bulk modulus K and higher G 0 lead to more than 15% difference from the baseline value. Increasing G 0 from 5.99e-3 MPa to 0.599 MPa led to a peak acceleration of 108g and a maximum Von Mises stress of 9.5 MPa compared to 79g and 12 MPa respectively with the original values. Results of this parametric study showed that brain constitutive law is found to have a minimal influence on skull deformation and stress as well as peak acceleration for the reconstructed drop test. Only skull mechanical parameters are observed in the parametric study since only the skull behaviour is investigated. However brain mechanical parameters should be studied if brain injuries are investigated with finite element models.

Application to a real-world head trauma

In order to illustrate the ability of the developed FE model to reproduce head injuries, a realworld head trauma was simulated with the new model. The accident was collected in the paediatric emergency department of a hospital. The selected case for accident reconstruction was a well-described accident case with precise explanation of the accident in the medical file. The description was made by the emergency medical service. Information required for numerical reconstruction were extracted from the medical files, e.g. fall height, nature of impact surface floor, area of impact and observed injuries. The accident involved a 1-month-old newborn child who fell on a concrete surface on the vertex. The head impacted first. A small bruise was observed on the vertex area of the head, which confirmed the area of impact. The origin of the fracture line and impact point was clearly identified by the medical staff and confirmed by medical imaging. For this accident, the impact point was identified in the vertex area leading to bilateral parietal fracture. The medical file clearly stated that the child did not lose consciousness. The observed injury was a skull fracture. Figure 5 illustrates fracture line resulting from this fall. CT images of the brain were used to conclude that no brain damage was observed. A medical consultation was performed 1 week after the fall, which indicated no neurological relapses. The nature of the impacted surface provides the mechanical characteristics of the floor.

Based on the height of the fall of 1 meter, which is stated in the medical file, an impact velocity was evaluated by kinetic and potential energy balance according to equation (1):

h g v . . 2  (1)
where h is the height of fall, and g is the gravity acceleration.

The FE model of the newborn was submitted to a direct impact against the modelled surface with an initial velocity of 4.4 m/s corresponding to 1m high of fall. The impacted surface was a concrete surface with a Young's modulus of 30 GPa and was modelled with shell elements with a thickness of 7mm.

Results in terms of skull fractures are given in Figure 6, showing a good accordance in terms of the fracture line. The fracture line of both parietal bones was observed with the model, and appears 5 ms after the impact. The length of the right and left parietal fracture was almost 67 mm and 85 mm respectively in the simulation.

It is also important to notice the dynamical behaviour of the skull at the time of impact which is an important characteristic of skull fracture in young infants, and quite different from adult fracture. The elasticity of the skull allows great malleability which lead to a "Ping Pong" skull response [START_REF] Mastrapa | Depressed skull fracture in Ping Pong: elevation with Medeva extractor[END_REF]. The impacted convex surface of the head becomes concave after a trauma. The numerical simulation of the accident reproduced this typical behaviour of the infant skull, as illustrated in Figure 7.

Discussion.

In this study, a biofidelic Finite Element Model of a newborn head was developed. Presently, the large number of existing FEMs of adult heads which are able to predict occurrence of neurological injuries, subdural haematomas or skull fractures, contrasts with the few existing paediatric head models available [START_REF] Desantis Klinich | Estimating infant head injury criteria and impact response using crash reconstruction and finite element modeling[END_REF][START_REF] Roth | Finite element analysis of impact and shaking inflicted to a child[END_REF][START_REF] Roth | Limitation of scaling methods in child head finite element modelling[END_REF][START_REF] Roth | Child head injury criteria investigation through numerical simulation of real world trauma[END_REF][START_REF] Coats | Parametric study of head impact in the infant[END_REF]. Experimental tests on cadavers provided data for validation of finite element models of the head. Concerning adult head, several experimental studies on brain pressure, displacement or skull fracture have been performed, and comparison between these tests and numerical results were performed in term of pressure (Nahum's et al [START_REF] Nahum | Intracranial pressure dynamics during head impact[END_REF]) or brain displacement (Hardy's et al. [START_REF] Hardy | Investigation of head injury mechanisms using neutral density technology and high-speed biplanar X-ray[END_REF]) and skull fracture (Yoganandan's et al. [START_REF] Yoganandan | Biomechanics of skull fracture[END_REF]) which ensured a realistic behaviour of the head when submitted to load cases. Moreover, these models also permitted the derivation of improved head injury criteria from an extensive real-world head trauma campaign (Deck et al. 2008 [29])

Validation

Concerning child head models, only Prange's experiments provided data on the global behaviour of the head, in terms of compressive force-displacement curve and acceleration of the head under impact. The proposed finite element model of the newborn head was validated against these experimental data. In the framework of the model validation, the global behaviour of the newborn skull was studied and its biofidelity was checked, as illustrated on Figure 3 and4. However, some discrepancies were observed in terms of accelerations curves. Several parameters could explain the difference of amplitude.

The newborn head undergoes wide and rapid changes during growth in terms of characteristic dimensions and mass. The present model represented a 50% newborn head as described in the Methods section. However, experimental tests performed by Prange et al. involved 1, 3 and 11-day-old specimens. This variation of ages of the tested cadavers could also lead to different skull responses in terms of acceleration since the mass of the head was slightly different.

Experimental tests have been performed with different cadavers having different mass (36% difference between specimens), which explain the large thickness of corridors. The model mass was almost 720 grammes compared to the experimental masses of 650 grammes. This difference of mass can be explained by the age of concerned heads (the FE model was based on a 17-day-old child whereas specimens are aged of 3 to 11-day-old) Furthermore, acceleration is usually evaluated at the centre of gravity (cog) of the head for rigid entities (like instrumented dummies' head). In this case, experimentation and simulation involved deformable structures, where no acceleration of the cog can be evaluated. Acceleration has been evaluated for both experimental tests and numerical simulations by dividing the impact force by the head mass. This way to calculate the head acceleration could be questioned and could lead to uncertainties concerning peak acceleration. Large standard deviations of peak acceleration in experimental tests illustrate this point.

Finally, skull fracture validation cannot be performed since the experimental tests do not investigate skull rupture. Reconstruction of well-documented real-world head trauma can help in the study of skull fracture. Even if information in the medical file do not allow a precise measurement of the length of fracture line, comparison between the simulation and the trauma can be made at a qualitative level: the simulated accident reproduces in a realistic manner the skull fracture line and the "Ping Pong" fracture, and confirms the biofidelity of the skull constitutive law.

Limitations and future work

Even if the model has shown a biofidelic behaviour, limitations concerning this study exist, especially concerning the material properties. Material properties of the brain tissue can be questioned. Even if brain properties are not the aim of this study, and even if its influence on skull response is very limited (parametric study), it is important to notice that only the linear viscoelastic aspect of the brain behaviour was taken into account, whereas more and more studies express its nonlinear aspect. Therefore, further investigations are needed to implement more biofidelic material laws for the brain tissue, especially if neurological injuries are studied.

Concerning skull properties, this study investigated skull deformation. The study did not provide skull fracture validation. Even if the accident reconstruction have shown interesting results in terms of rupture, fracture was only shown by the real-world trauma, and not at an experimental level, as it was performed on adult head [START_REF] Yoganandan | Biomechanics of skull fracture[END_REF]. In that case, precise experimental tests on skull fracture would be necessary, and it is obvious that a great number of similar accident cases are needed and will be investigated in order to consolidate the rupture behaviour of the skull.

Conclusion

A new finite element model of the newborn head was developed in this study, and was validated against experimental data in terms of skull deformation and acceleration. Child biomechanics suffers strongly from great limitations in terms of experimental data. Ethical issues limit the type of experiments which can be conducted, yet these experiments are required for mathematical model validation. To the author knowledge, only Prange's test allowed an acceptable global validation of infant head impact under dynamic compression and impact. These data were used to validate the new FE model in terms of skull deformation and acceleration. The developed model has been used to reconstruct a real-world head trauma leading to encouraging results in terms of skull fracture and "ping pong" skull deformation. Previously developed child head FE models have been used to investigate paediatric head trauma, but none of these has been validated due to the lack of experimental data for validation. As for adult's head, validation against experimental data is essential in the development of efficient injury predicting tools, but remains one of the main difficulties in FE studies concerning the paediatric population. With this first step in terms of FE validation, this study presents a child FE head model which, with further development, could be used for analysing paediatric head trauma in future