Notations and conventions

• We will work in 2D.

• upercase letters will always be vectors.

• lowercase are scalars or scalar functions • . • . or < ., . > unconsistently represent either scalar product of vectors or matrix vector multiplications. My lousy mathematical education is to blame for this but I really like using both ... and in different context ... I find the result more readable ... I am sorry.

• Numerical results are best viewed in color.

Introduction

The Eikonal equation and its companion the Transport equation arise from WKB (high frequency) asymptotics of various Linear and Non-linear wave equations (see [START_REF] Fedoryuk | Partial Differential Equations[END_REF] for a nice review of WKB methods). As a natural front propagator, it also has a strong relationship with the level set methods with many applications, in particular in image processing ( [START_REF] Tsai | Total variation and level set methods in image science[END_REF] [3]...) I think this probably is the shortest motivation section I ever wrote ...

Locally, the Eikonal equation simply prescribes the gradient of the unknown scalar function φ to be on the sphere of radius n(X) at each point X. ∇φ(X) = n(X).

(1)

Obviously, this equation makes sense for positive right hand side and local classical solutions only exist for strictly positive smooth (C 1 ∪ W 2,∞ ) given functions n, an assumption we will stick to in this paper. Assuming that the gradient (or a gradient) exists and has at least Lipschitz derivatives on some open set, one can construct a solution by the method of characteristics : let us define the integral curves of the gradient Y (s) (also called "rays")

d ds Y (s) = ∇φ(Y (s)) (2)
and set P = ∇φ(Y ). We easily eliminate φ using (1)

dP ds = D 2 φ • d ds Y = D 2 φ • ∇φ(Y ) = 1 2 ∇n 2 (Y (s)) (3)
The ODE system (2) (3) can be solved independently from [START_REF] Fedoryuk | Partial Differential Equations[END_REF] and the "phase" function φ can be computed along the Lagrangian coordinate Y using a third ODE

d ds φ(Y (s)) = dY ds • ∇φ(Y (s)) = P (s) 2 = n 2 (Y (s)) (4)
It gets more complicated when it comes to global solutions, for instance : solving (1) in some domain Ω or Ω c (exterior problem, Ω can be just a source point) with Dirichlet boundary conditions φ = φ 0 , on ∂Ω.

(

) 5
The existence of a unique smooth gradient field compatible with the boundary conditions (even if there are themselves also smooth) is not garanteed. Even worse : it is fairly easy to construct many weak solutions by patching local solutions obtained by the method of characteristics. In this paper we will focus on the so called viscosity solution of this equation wich ensures that (1)( 5) is well posed. The viscosity solution φ v can conveniently be characterized as minimization over all paths Y of the optical length (this is in essence Fermat principle) :

φ v (X) = inf 8 > > < > > : Y 0 , Y ∈ W 1,∞ (0, t) Y (t) = X Y (0) = Y 0 φ 0 (Y 0 ) + t 0 n(Y (s)) Ẏ (s) ds (6)
Mathematically the solutions of (2)(3) (with ad-hoc initial conditions) are the Euler Lagrange equations of [START_REF] Abgrall | Construction of simple, stable and convergent high order schemes for steady first order Hamilton Jacobi equations INRIA tech[END_REF]. The viscosity solution is therefore the minimum "phase" value out of all possible solutions given by ( 4). This produces bounded uniformly continuous solutions that are generally C1 except possibly along curves (or surfaces on 3D) called kinks where the gradient is discontinuous. This minimum property contributes to the well posedness of [START_REF] Fedoryuk | Partial Differential Equations[END_REF] in this class of solution and also helps making upwind finite difference aproximations convergent. This is further discussed in section 2 and 3 where the widely used Rouy Tourin (RT) Finite Difference Scheme [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF] is reviewed and analysed in terms of Lagrangian characteristics.

As made explicit in the computations above, the method of characteristics, a classics for linear hyperbolic problems, can only be applied to the non-linear Eikonal equation (1) because of its Lagrangian structure :

1. As already mentionned [START_REF] Fedoryuk | Partial Differential Equations[END_REF] constrains the gradient to be on a sphere and solving locally the Eikonal equation reduces to find an angle.

2. Assuming the second order derivatives exists, (1) also yields

D 2 φ • ∇φ = 1 2 ∇n 2 (7)
which express that the derivatives of n prescribe the second derivatives of φ in the Lagrangian trajectory direction.

This paper is a discussion on the impact of this Lagrangian structure in the understanding and the design of Eulerian numerical methods for the Eikonal equation. In particular we uncover a superconvergence phenomena of RT algorithm which was not apparently known

INRIA

and also derive an elegant compact upwind second order scheme.

Readers interested in the theory of viscosity solutions can find more details in [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] [5] Let us quickly mention that this theory was originally develloped for a general class of first order Hamilton-Jacobi equations

H(X, ∇φ(X)) = 0 ( 8
)
where the Hamiltonian function H(X, P ) is convex (or concave) and coercive in P , thus retaining the essential features the Lagrangian structure. in particular [START_REF] Osher | High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF] becomes

H x (X, ∇φ) + D 2 φ • H p (X, ∇φ) = 0 (9)
where H p is strictly monotonic and thus invertible. We have not written all the details but we believe that most of the discussion carried out in this paper apply to (8).

The Rouy-Tourin (RT) Algorithm

We briefly review the RT first order Finite Difference scheme [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF] from a Lagrangian view point.

The first step consist in a rewrite of (1) as a pointwise optimization problem :

∇φ(X) = n(X) ⇔ sup Q ≤1 {∇φ(X) • Q -n(X)} = 0. (10)
This is easily proven by writing both ∇φ and Q in polar coordinates. The optimal Q is given by Q opt = ∇φ ∇φ .

Then, discretize the directional derivative by introducing a fictitious "time" :

sup Q ≤1 { φ(X -dt Q) -φ(X) -dt -n(X)} = 0 ⇔ φ(X) = inf Q ≤1 φ(X -dt Q) + dt n(X) (11
) Remark that we have slightly manipulated the expression to make it appear consistent with [START_REF] Abgrall | Construction of simple, stable and convergent high order schemes for steady first order Hamilton Jacobi equations INRIA tech[END_REF]. Indeed remember that 1. Q opt = ∇φ ∇φ is the direction of a ray passing through X, see (2).

2. dt n(X) is a local approximation of the optical length.

Only by minimizing the phase over locally over "portions of rays" can we recover the minimum over global rays initialized on the boundary, i.e. the viscosity solution. Conversely only an absolute minimum can be computed by local minimization steps. This unfortunately makes plain viscosity solutions irrelevant for the multivalued traveltime problem, for more

X Ray direction 0,!h 0,h 0,0 !h,!h !h,0 h,0 Upwind quadrant Y Figure 1:
Stencil and upwind quadrant on this topic see [11] [9]).

In the next section we show that RT algorithm can indeed be linked to ray tracing but let us first finish the derivation of the algorithm. We follow a slightly different technique which we find simpler than the presentation in [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF]. Let h be a discretization step in x and y the cartesian coordinates X = {x, y}, ({0, 0}, {h, 0}, {0, h}, {-h, 0}, {0, -h}) are the coordinates of the usual regular grid five point stencil around the origin (figure 1) and (φ 0,0 , φ h,0 , φ 0,h , φ -h,0 , φ 0,-h ) are the unknown values of the phase at these grid points. We first assume that the direction of the "optimal" ray passing through {0, 0} points in the x ≥ 0, y ≥ 0 quadrant (figure 1) . This means that the point X -dt Q opt should be chosen in the x ≤ 0, y ≤ 0 quadrant, this is know as upwinding.

We focus on one this quadrant, denote ∇φ = (a, b) and Taylor expand

φ 0,-h = φ 0,0 -ha + O(h 2 ) φ -h,0 = φ 0,0 -hb + O(h 2 )
As usual we make a linear approximation by neglecting the second order error terms, the Eikonal equation becomes

a 2 + b 2 = ( φ 0,0 -φ 0,-h h ) 2 + ( φ 0,0 -φ -h,0 h ) 2 = n 0,0

INRIA

Note that the assumption that the ray points in the x ≥ 0, y ≥ 0 quadrant is equivalent to φ 0,0 -φ 0,-h ≥ 0 and φ 0,0 -φ -h,0 ≥ 0 : the phase increases in the direction of the ray. This remark combined with the using the inf property of (11) leads to the four quadrant formula

max(a + l , a - r ) 2 + max(b + l , b - r ) 2 = n 0,0 (12)
where a + = max(0, a), a -= max(0, -a) and

a l = φ 0,0 -φ -h,0 h a r = φ 0,0 -φ h,0 h b l = φ 0,0 -φ 0,-h h b r = φ 0,0 -φ 0,h h
We end the discretization of (1)(5) by writing ( 12) at all interior grid points and using the Dirichlet boundary conditions. This yields a system of non-linear equations for (φ i,j ) i,j , the phase at all grid points. This can more abstractly written as g(φ i,j , φ i+1,j , φ i-1,j , φ i,j+1 , φ i,j-1 ) = 0 [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF] for all interior grid points (i, j), boundary point values are specified using the Dirichlet boundary condition. The function g is the called the numerical Hamiltonian and embodies formula [START_REF] Qian | An adaptive finite difference method for travel-time and amplitudes[END_REF]. Existence and uniqueness of the discrete solution is proved in [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF]. It theoretically is O( √ h) accurate but it is also well know that the expected first order accuracy is the rule away from the kinks where the gradient is discontinous. The impact of the regularity of the boundary conditions (especially point sources) also needs attention (see [START_REF] Qian | An adaptive finite difference method for travel-time and amplitudes[END_REF]).

We finally adress the problem of solving this non-linear system. Once again the Lagrangian interpretation of scheme [START_REF] Engquist | Multiphase computation in geometrical optics[END_REF] is helpful. Information clearly propagate along the Lagrangian wind and the value at the center of the stencil is "down the local wind". It therefore makes sense to relax [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF] and compute iteratively a sequence of solutions (φ k i,j ) i,j as g(φ k+1 i,j , φ k i+1,j , φ k i-1,j , φ k i,j+1 , φ k i,j-1 ) = 0 [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] for all interior grid points (i, j). Getting back to our simple one quadrant description this is easily solved using

φ k+1 0,0 = 0.5 (φ k 0,-h + φ k -h,0 ) + 0.5 2 h 2 n 2 0,0 -(φ k 0,-h -φ k -h,0 ) 2 Assuming (φ k 0,-h -φ k -h,0 ) 2 < h 2 n 2 0,0 φ 0,0 = min(φ k -h,0 , φ k 0,-h ) + h n 0,0 else
Once again we can generalized to the four quadrant using the inf property

α = min(φ k 0,-h , φ k 0,h ) β = min(φ k -h,0 , φ k h,0 ) φ k+1 0,0 = 0.5 (α + β) + 0.5 2 h 2 n 2 0,0 -(α -β) 2 Assuming (α -β) 2 < h 2 n 2 0,0 φ 0,0 = min(α, β) + h n 0,0 else (15)

RR n°0123456789

The first minimization pick up the upwind quadrant, the last test is then just activated when the optimal direction is exactly on a grid line. This iterative process has been shown to produce a monotonicaly convergent sequence to the solution of [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF]. The grid point ordering of the update ( 14) for each iteration k has been much studied and has a very important impact on the speed of convergence, i.e. the number of iterations needed for [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] to converge to [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF]. An "optimal" strategy clearly needs to take into account some kind of Lagrangian information. The famous fast marching method [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF] achieves a O(N logN ) complexity , N beeing the number of grid points. at the expense of a sophisticated sorting algorithm. Group marching algorithm [START_REF] Kim | Seongjai An O(N ) level set method for eikonal equations[END_REF] is elegant even though it does not seem to be as simple and efficient of the Sweeping algorithm [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF] [START_REF] Kim | Seongjai An O(N ) level set method for eikonal equations[END_REF] which is in essence just Gauss-Seidel combined with the inf property and alternate choice of "sweeping" direction {x > 0, x < 0, y > 0, y < 0}. Sweeping seems in most cases to be of complexity O(N ) and was proven to be so for n ≡ 1 [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]. We rediscuss sweeping as a motivation in the next section and also in section 6 where we propose our own implementation of the Sweeping method.

Eulerian based Ray tracing (EBRT)

As explained in the introduction, the method of characteristics (ray tracing) can be considered as the fundamental method of resolution of the Eikonal equation. It also plays a role through upwinding (see 11) in designing Eulerian numerical methods such as the RT algorithm. In this section we try to walk our way backward from the discrete viscosity solution (φ i,j ) i,j of the preceding section to the rays underlying the solution, that is the optimal paths of problem [START_REF] Abgrall | Construction of simple, stable and convergent high order schemes for steady first order Hamilton Jacobi equations INRIA tech[END_REF]. The starting motivation was, other than mere curiosity, linked to the Sweeping method for which there are no rigorous estimate on the number sweep needed to converge. One possible understanding of sweeping method is that it picks portions of correct rays along the direction is propagates and that it therefore takes as many iterations of sweeps in all directions than 45 o turns in the underlying exact rays to propagate the correct phase from portions to portions. So it seems that a reasonable first step towards a convergence speed theory is to make sure that we can really interpret locally the discrete solution as local ray tracing : we want to check that ( 11) is really a local discrete version of ( 6). So we propose to construct "discrete rays" from RT discrete Eulerian solution by the following EBRT procedure :

1. Compute the viscosity solution φ i,j on a regular grid.

2. Build φ h a continuous P1 linear interpolant of φ i,j .

3. From any point X do a gradient descent on the surface. Note that step 2. is not uniquely determined as for each group of four neighbors points we can split the cell quadrangle along either diagonal to build φ h . Numerical experiment do not seems to be sensitive to this issue which however raises a technical difficulty (discussed 12) the remaining point is attached to another neighboring cell. In order to make the discussion as simple as possible we choose the diagonal which is consistent with the scheme. Steps 1. to 3. produce Y h (s) a polygonal path solution of Ẏh = -∇φ h for all (τ i ) i∈I , I being the set of indices of triangles crossed by Y h . Because there are no sinks (φ i,j cannot be strictly smaller than all its neighbors) and no plateau in the φ h surface (n h is strictly positive), Y h must exit the domain in a finite number of triangles, say at Y 0 . There is however no garantee of uniqueness for Y h in particular if it hits a grid point where branching could occur. We resort once again to the inf property and choose the steepest descent direction in that case. Figure 2 and3 show this algorithm in action both for n ≡ 1 (homogenous index of refraction) and n as a Gaussian, thus bending the rays (heterogeneous). A zoom details the φ h triangles used. We also plotted the contour lines of the phase.

INRIA
We now want to check that as h goes to 0 our EBRT solution converges to the actual piece of ray yielding the viscosity solution. Let's focus on one triangle τ i , φ h is built to satisfy ∇φ h = n h where n h is the value at the center grid point of the scheme (this is not true on all triangles as discussed below). Let Y h : Y in → Y out be the piece of path going through our triangle The phase difference satisfy

δφ i = φ h (Y in ) -φ h (Y out ) = n h Y out -Y in = Yin Yout n h (Y h (s))ds (16)
and as n h is a piecewise constant approximation of n

δφ i = Yin Yout n(Y h (s)) + O(h)ds = Yin Yout n(Y h (s))ds + O(h 2 ) RR n°0123456789
φ h (X) = φ 0 (Y 0 ) + i∈I δφ i = φ 0 (Y 0 ) + X Y0 n(Y h (s))ds + O(h)
We now recall that φ h → φ v uniformly (standard convergence results are not established for the P 1 interpolant but φ h can be bounded by convergent sub and super solutions) and

(Y h ) h is uniformly bounded in C 0 ∪ W 1,∞ . So, for any convergent subsequence Y h ′ → Y ′ : φ v (X) = φ 0 (Y 0 ) + X Y0 n(Y ′ (s))ds
We conclude using the inf property in [START_REF] Abgrall | Construction of simple, stable and convergent high order schemes for steady first order Hamilton Jacobi equations INRIA tech[END_REF]. The optimal path is necessarily a ray, uniquely determined by intial and end points.

There is a technical issue linked to step 2. of our EBRT procedure : ∇φ h = n h may not be satisfied exactly on all (τ i )s.

We then distinguish two cases :

1. ∇φ h < n h . Then, it is fine because

δφ i = φ h (Y in ) -φ h (Y out ) < n h Y out -Y in = Yin Yout n h (Y h (s))ds 2. ∇φ h > n h . Then, instead of gradient descent choose Y out as Y out = Y in + e θ Y out -
Y in note that Y out depends on the angle θ. Then

δφ i = φ h (Y in ) -φ h (Y out ) = ( e θ • ∇φ h ∇φ h ) ∇φ h Y out -Y in

INRIA

If ∇φ h is not too far from n h (looks likely) then we can pick θ such that ( e θ • ∇φ h ∇φ h ) ∇φ h = n h . There is two possible angle on each side. Use again inf principle to make this step uniquely determined and we recover [START_REF] Kim | Seongjai An O(N ) level set method for eikonal equations[END_REF].

Summing again and passing to the limit we get this time

φ v (X) ≥ φ h (Y 0 ) + X Y0 n(Y ′ (s))ds
and can conclude exactly as before as the inf property ( 6) actually gives us equality in the above relation.

We now check the convergence numerically. We give in the tables below both the phase error at the starting point (line Phi) the error on the exit point between the exact ray and the EBRT (line EBRT) output as a function of N the number of points used along one dimension (h = Both tables show the expected first order convergence on the phase function and a similar convergence for the EBRT procedure.

Superconvergence of upwind schemes

Section 3 leaves us with a paradox. If RT algorithm is first order and our EBRT procedures consists in solving d ds Y (s) = ∇φ h (Y (s)). How is possible to get convergence ? Indeed ∇φ h should be O(1) with respect to h and not convergent.

It is worth noticing that the gradient does not really appear in the convergence arguments of section 3. An indeed a simple numerical study corroborates that the gradient is first order : We computed the max error for a smooth solution (homogeneous space experiment figure 2) both for the phase (line Phi) and the x derivative of the phase (line Phi x). The table below indicated these values for different discretization N (h = First order convergence (O(h)) can be observed both for the phase and its derivatives.

We make two remarks we think are relevant to understand this superconvergence phenomena.

First we recall that we made a first order approximation in [START_REF] Engquist | Multiphase computation in geometrical optics[END_REF] sup

Q ≤1 { φ(X -dt Q) -φ(X) -dt -n(X)} = O(dt 2 )
If we instead use a second order Taylor expansion we get

sup Q ≤1 { φ(X -dt Q) -φ(X) -dt + dt 2 < Q, D 2 φ(X) • Q > -n(X)} = O(dt 3 ) Remembering Q opt = ∇φ ∇φ = ∇φ n . We obtain    φ(X -dt ∇φ(X) n(X) ) -φ(X) -dt + dt 2 < ∇φ(X) n(X , D 2 φ(X) • ∇φ(X) n(X) >    -n(X) = O(dt 3 ) But (7) : D 2 φ(X) • ∇φ(X) = 1 2 ∇(n 2 (X)) = n(X)∇n(X)
. So we can rearrange to recover a formula similar to [START_REF] Engquist | Multiphase computation in geometrical optics[END_REF] where a trapezoidal rule is used to approximate the local optical length :

φ(X) = φ(X -dt ∇φ(X) n(X) ) -dt dt 2 < ∇φ(X) n(X) , ∇n(X) > -n(X) + O(dt 2 ) φ(X -dt ∇φ(X) n(X) ) -dt 1 2 < dt ∇φ(X) n(X) , ∇n(X) > -n(X) + O(dt 2 ) φ(X -dt ∇φ(X) n(X) ) + dt 1 2 n(X) + n(X -dt ∇φ(X) n(X) ) + O(dt 2 )

INRIA

Apparently second order upwinding only requires a better approximation n and this is in particular alway true for constant indices of refraction. Of course, in deriving RT, we still have the piecewise linear approximation of φ(X -dt Q opt ) which produces a first order error.

A second Order Scheme

The most popular high order schemes for Hamilton-Jacobi equation are based on the ENO (Essentially Non Oscillatory) technique [START_REF] Osher | High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF] which automatically moves the second order stencil around a grid point to fetch the smoother (kink free) approximation. The approach followed here is different. Using the Lagrangian structure of the Eikonal equation we propose a pure upwind second order extension of R-T scheme that needs only one more point added to the stencil.

We first write a second order Taylor expansion

φ(X + δX) = φ(X) + δX • ∇φ(X) + 1 2 < δX, D 2 φ(X) • δX > +O(δX 3 ) (17)
We want to take advantage of ( 7) : D 2 φ(X) • ∇φ(X) = n(X)∇n(X) and decompose the variation δ in the local coordinate system made of the Lagrangian direction and its orthog-

onal δX = 1 n 2 (X) (δX • ∇φ) ∇φ + (δX • ∇φ ⊥ ) ∇φ ⊥
Plugin this decomposition in the second order term of (17) yields

1 2 < δX, D 2 φ(X) • δX >= 1 2 n 4 (δX • ∇φ) 2 < ∇φ, D 2 φ • ∇φ > + 1 n 4 (δX • ∇φ)(δX • ∇φ ⊥ ) < ∇φ ⊥ , D 2 φ • ∇φ > + 1 2 n 4 (δX • ∇φ ⊥ ) 2 < ∇φ ⊥ , D 2 φ • ∇φ ⊥ > (18)
We first restrict to the simpler homogeneous case n ≡ 1. Derivatives of n vanish and (17) becomes

φ(X + δX) = φ(X) + δX • ∇φ(X) + 1 2 (δX • ∇φ ⊥ ) 2 < ∇φ ⊥ , D 2 φ • ∇φ ⊥ > +O(δX 3 ) (19)
For point sources solutions φ(X) = X -X 0 , < ∇φ ⊥ , D 2 φ • ∇φ ⊥ > is the curvature of the level sets of the phase (also called fronts). We will denote it C. Our new approximation of the phase (19) uses "shape function" of the form (denote ∇φ = (a, b))

φ(X + δX) = φ(X) + δX • (a, b) + C 2 n 4 (δX • (-b, a)) 2 + O(δX 3 ) (20) RR n°0123456789
This is a parabolic "local curvature" correction to RT piecewise linear interpolation (section 2).

We now proceed as in the RT scheme but using the second order approximation (20). We first restrict to a quadrant. To realize our approximation we need to determine the three coefficients (a, b, C). We therefore need three equation and we use the four grid points ({0, 0}, {-h, 0}, {0, -h}, {-h, -h}) to obtain

           φ -h,0 = φ 0,0 -h a + C 2 (-h b) 2 + O(h 3 ) φ 0,-h = φ 0,0 -h b + C 2 (h a) 2 + O(h 3 ) φ -h,-h = φ 0,0 -h (a + b) + C 2 (h (b -a)) 2 + O(h 3 ) (21)
Truncating the third order error we can get a grid value approximation of the curvature C and eliminate the last equation for instance.

line 1 + line 2 -line 3 ⇔ C = φ 0,-h -φ 0,0 + φ -h,0 -φ -h,-h h 2 a b (22)
The next natural steps should be 1. to solve for (a, b) as a function of (φ 0,0 , φ -h,0 , φ 0,-h , φ -h,-h ) and replace in the Eikonal equation to write the numerical Hamiltonian.

2. relax : compute φ 0,0 as a function of the values at grid points ({-h, 0}, {0, -h}, {-h, -h}).

This turn our to be a difficult non linear problem, involving the roots of a sixth order polynomial. Before explaining how to simplify this into a more tractable problem. We derive a natural and elegant interpretation of our scheme. Remembering the Eikonal equation can be recast locally as {a, b} = {n cos(θ), n sin(θ)}, we can use (21)(22) to write the a local weighted centered FD formula in the upwind quadrant

   a = (1 -1 2 tan(θ)) φ 0,0 -φ -h,0 h + 1 2 tan(θ) φ 0,-h -φ -h,-h h b = (1 -1 2 tan -1 (θ)) φ 0,0 -φ 0,-h h + 1 2 tan -1 (θ) φ -h,0 -φ -h,-h h (23)
This is not a practical solution formula as θ depends on (a, b) but it shows that we average the local first order finite difference (FD) according to the Lagrangian direction. For instance if the ray makes a π 4 angle with the x axis, i.e. exactly cuts in half the quadrant, then we use an arithmetic average of the FD derivatives of the four grid points of this quadrant. If the ray approaches an axis then the weight are non-linearly adjusted by the tan(θ) functions.

In order to simplify our local solver we remember that RT scheme 1. determines the upwind quadrant.

INRIA

2. does so with first order accuracy.

We will denote the first order gradient obtained from RT as (a 1 , b 1 ) and assume the gradient points in the positive x and y direction. We can therefore work with the upwind quadrant ({0, 0}, {-h, 0}, {0, -h}, {-h, -h}). This also means that

a1 = φ 1 0,0 -φ -h,0 h b1 = φ 1 0,0 -φ 0,-h h ( 24
)
are given as a function of φ 1 0,0 the first order RT solution. Using section 4 remark that (a 1 , b 1 ) = (a, b) + O(h), we notice that replacing in the quadratic term of the Taylor expansions (21) preserves second order accuracy

             φ -h,0 = φ 0,0 -h a + C 2 (-h b 1 ) 2 + O(h 3 ) φ 0,-h = φ 0,0 -h b + C 2 (h a 1 ) 2 + O(h 3 ) φ -h,-h = φ 0,0 -h (a + b) + C 2 (h (b 1 -a 1 )) 2 + O(h 3 ) (25)
We can again eliminate C with formula ( 22) where (a, b) is replaced by (a 1 , b 1 ). Then, (23) simplifies to

   a = (1 -1 2 tan(θ 1 )) φ 0,0 -φ -h,0 h + 1 2 tan(θ 1 ) φ 0,-h -φ -h,-h h b = (1 -1 2 tan -1 (θ 1 )) φ 0,0 -φ 0,-h h + 1 2 tan -1 (θ 1 ) φ -h,0 -φ -h,-h h (26)
where θ 1 is the first order angle given by RT and therefore tan(θ 1 ) = b1 a1 . We can now write our second order numerical Hamiltonian (restricted to the upwind quadrant)

g(φ 0,0 , φ 0,-h , φ -h,0 , φ -h,-h ) = a 2 + b 2 -n 0,0 = 0 (27)
where a, b is given by (26).

About the properties of this numerical scheme :

1. Consistency of order 2.

2. Stability as long as θ 1 stays away from grid lines.

3. Monotonicity, as a second order scheme is of course non monotone (it is easily checked that g can not be simultaneously a decreasing function of φ 0,-h and φ -h,0 at the same time).

So the usual convergence theory [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] does not apply here. However a closer look (in the spirit of [START_REF] Abgrall | Construction of simple, stable and convergent high order schemes for steady first order Hamilton Jacobi equations INRIA tech[END_REF] th. 3.1 1 to the classic viscosity convergence shows that it is enough to be a O(h) perturbation of a monotonic scheme to be convergent. We notice that (26) can also be written

   a = a 1 -1 2 tan(θ 1 ) { φ 0,0 -φ -h,0 h - φ 0,-h -φ -h,-h h } b = b 1 -1 2 tan -1 (θ 1 ) { φ 0,0 -φ 0,-h h - φ -h,0 -φ -h,-h h } (31)
so as long as the tan(θ 1 ) remains bounded above and from 0, and independently of h (which means that the ray direction should be away from an axis) any smooth grid sampled C 2 function φ will satisfy

g(φ 0,0 , φ 0,-h , φ -h,0 , φ -h,-h ) = a 2 1 + b 2 1 -n 0,0 + O(h) ( 32
)
and it is well known that a 2 1 + b 2 1 leads up to a monotonic numerical Hamiltonian.

For the upper right quadrant, a condition like (see figure 4)

2 -ǫ > θ 1 > ǫ (33)
where ǫ > 0 is independent of h is sufficient to garrantee convergence of the second order scheme.

Of course, there is no way to systematically enforce (33) as it depends on the problem data. We must keep away from the second order scheme when (33) is not satisfied. In that case (assume for instance that θ 1 < ǫ) b 1 = sin(θ 1 ) = O(ǫ), a 1 = 1 -O(ǫ) , the upwind Taylor expansion along the x axis gives

φ -h,0 = φ 0,0 -h a + h O(ǫ) + C 2 (-h b 1 ) 2 + O(h 3 ) = φ 0,0 -h + h O(ǫ) + O(h 2 ) O(ǫ) 2 + O(h 3 ) ( 34
)
1 In the classic proof the key point is to consider sequences of approximations hn → 0 , Xn → x 0 , φ hn (Xn) → φ(X 0 ) and show that φ is a viscosity super-solution of the Eikonal equation. To do so one assumes that X 0 is a point of local minimum for φψ and the above sequences are built such in such a way that Xn is a global minimum of φ hn -ψ in a neighborood of X 0 independent of h. Then one sets ξn = φ hn (Xn) -ψ(Xn) and gets ξn → 0 , φ hn ≥ ψ + ξn near X 0 . This is where monotonicity of the numerical Hamiltonian comes in. Replacing the φ hn values by ψ + ξn we get (see [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF] H num (ψ + ξn) = g(ψ i,j + ξn, ψ i,j+1 + ξn, ψ i+1,j + ξn, ψ i-1,j + ξn, ψ i,j-1 + ξn) ≥ 0 (28)

Then pass to the limit (remark that in this first order equation the ξn cancel immediately) to get, using consistency (H is the Eikonal equation Hamiltonian).

H(X, ∇ψ(X)) ≥ 0 (29) which proves the viscosity super-solution of φ.

In [START_REF] Abgrall | Construction of simple, stable and convergent high order schemes for steady first order Hamilton Jacobi equations INRIA tech[END_REF] th. 3.1, the author uses a clever blending of first order monotonic and high order numerical Hamiltonians to satisfy approximately (28) and still get (29) at the limit. In our case this is much simpler as we get H num (ψ + ξn) + O(hn) ≥ 0 (30) where H num is Rouy Tourin first order monotonic Hamiltonian, and can easily pass to the limit (29). If we use this formula to update φ 0,0 we get hO(ǫ) accuracy. In practice we just use the already computed first order value φ 1 0,0 .

INRIA ! Y X " # ! # 2 /
Remark : Of course this treatment will only hold as long as ǫ = O(h) and ǫ need to fixed independent of h p.15 16 on the convergence of the scheme). So we will loose second order accuracy when h gets below a threshold depending on ǫ. The most reasonable treatment when θ 1 < ǫ is to change the stencil, the ray direction is along the x axis. To get a better approximation of the curvature of the front we can add the point (0, -2h) instead of (-h, -h). This will change formulae (26) and eliminate the problematic division by cos(θ 1 ) in the coefficients. This is currently under investigation ... [10].

Before summarizing the complete scheme we need to adress the resolution of the discrete system. As for the RT algorithm we use relaxation. We need to be able to solve for φ k+1 0,0

g(φ k+1 0,0 , φ k 0,-h , φ k -h,0 , φ k -h,-h ) = 0 (35)
This boils down to finding the roots of a second order polynomial. For simplicity we omit the k upperscripts, the equation is

a 2 + b 2 = n 2 0,0
where (a, b) are given in (31). This can be rewritten (carefully) as

A φ 2 0,0 + B φ 0,0 + C = 0 (36)
where (A, B, C) are easily computed and depend on the values at the other grid points and on (a 1 , b 1 ). When (33) is satisified, we expect (it works numerically and I need to find the energy to do the computation) (36) to have at least a real root close to φ 1 0,0 .

We recapitulate the algorithm in the constant index of refaction case RR n°0123456789

1. Iterate the relaxed RT first order Algorithm (14).

2. Make the second order correction (36) in the relevant quadrant if (33) is satisfied.

3. Otherwise, stay with the first order solution (34) else.

The heterogeneous index of refraction works the same but the algebra is more complicated. We keep all the terms in (18) but use [START_REF] Osher | High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF] to simplify, equation (19) becomes

φ(X + δX) = φ(X) + δX • ∇φ(X) + C 2 n 4 (δX • ∇φ ⊥ ) 2 + 1 n 3 (δX • ∇φ)(δX • ∇φ ⊥ )(∇φ ⊥ • ∇n) + 1 2 n 3 (δX • ∇φ) 2 (∇φ • ∇n) +O(δX 3 ) (37)
and ( 21), with the same notations augmented with n = n 0,0 , n x = ∂ x n 0,0 , n y = ∂ y n 0,0 and ∇φ ⊥ = (-b, a) corresponds to

               φ -h,0 = φ 0,0 -h a + h 2 2 n 4 (C b 2 + n a 2 (a n x + b n y ) -2 a b (-b n x + a n y )) φ 0,-h = φ 0,0 -h b + h 2 2 n 4 (C a 2 + n b 2 (a n x + b n y ) + 2 a b (-b n x + a n y )) φ -h,-h = φ 0,0 -h (a + b) + h 2 2 n 4 (C (a -b) 2 + n (a + b) 2 (a n x + b n y ) + 2 (a 2 -b 2 ) (-b n x + a n y )) (38
) Where we've truncated to third order accuracy. Again we can get a grid value approximation of the curvature C and eliminate the last equation for instance.

C = 1 a b ( n 4 h 2 (φ 0,-h -φ 0,0 + φ -h,0 -φ -h,-h ) + n n x a 3 + n n y b 3 ) (39)
We can then eliminate C and the last equation, set

F = φ 0,-h + φ -h,0 -φ -h,-h ...        φ -h,0 = (1 -b 2 a ) φ 0,0 -h a + b 2 a F + h 2 2 n 4 n n x (a 3 + 2 a b 2 + b 4 a )) φ 0,-h = (1 -a 2 b ) φ 0,0 -h b + a 2 b F + h 2 2 n 4 n n y (b 3 + 2 b a 2 + a 4 b )) (40)
Note that it can be simplified further into the remarkably simple

     φ -h,0 = (1 -b 2 a ) φ 0,0 -h a + b 2 a F + h 2 2 a n n x φ 0,-h = (1 -a 2 b ) φ 0,0 -h b + a 2 b F + h 2 2 b n n y (41)

INRIA

We have not (yet ?) been able to use this formulation so we go back to (40) and as in the homogeneous case this system with the additional Eikonal equation ({a, b} = n {cos(θ), sin(θ)}) must be solved for (a, b, φ 0,0 ). It is hopelessly non-linear. We thus propose as in the homogenous case to get back to a linear system in (a, b) by reverting to the first order quantities (a 1 , b 1 ) in the high order terms

       φ -h,0 = (1 -b 1 2 a 1 ) φ 0,0 -h a + b 1 2 a 1 F + h 2 2 n 4 n n x (a a 2 1 + 2 a 1 b 1 b + b b 3 1 a 1 )) φ 0,-h = (1 -a 1 2 b 1 ) φ 0,0 -h b + a 1 2 b 1 F + h 2 2 n 4 n n y (b b 2 1 + 2 b 1 a 1 a + a a 3 1 b 1 )) (42)
The first remark is that again we obtain a first order perturbation R-T numerical Hamilas these formula only differ from (26) by a O(h 2 ) perturbation. After relaxation the problem reduces to finding the roots of a second order polynomial. We pick up the root closest to the solution given by the first order scheme. We leave the courageous reader to the remainder of the algebra and turn to the numerical study.

Our first test case are cooked up to ensure that rays will propagate along the x = y axis and never have angles over ± π 4 with respect to this principal direction. The use of a staggered grid allows to consider always the same upwind quadrant (but one must be damn careful with the coordinate system in the implementation) see figure 5. Initial (boundary conditions) must be given on a 2 grid points layer.

Just one sweep of the grid in the direction of the rays is sufficient to compute the discrete solution. The table below show the max error at final depth for a constant index of refraction and a curved front (x10 It indeed confirms second order convergence.

In the next experiment the index is formed of a convex/concave lense and given by ((x, y) ∈ [-1, 1] 2 ) n(x, y) = 2.8 if x 2 + y 2 ≥ 0.5 n(x, y) = 2.8 ± 0.5 exp -30 (x 2 +y 2 ) else. In figures 6 7 we plot the phase at final depth for 1. RT first order scheme using 40 (black solid line) and 1600 (red solid line) grid points (on each axis). The red line can be considered as the exact solution.

Main axis of propagation X Y The second order convergence is obvious. Note how the kink is well resolved for the convex lense even with the second order scheme which stays purely upwind.

More numerical results the next section.

6 Notes on a Sweeping implementation (staggered grid)

In order to be able to take into account any ray direction we alternate sweeps on our staggered grid as indicated in figure 8.

Sweeping for first order scheme usually update the value given by the current sweep only if it gives a smaller value than the previous sweep (the phase is initialized to a very large value except at sources and boundary conditions). This is an efficient want to iteratively take into account the minimization over the four quadrant in [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]. Remark also that if the sweep determines the upwind quadrant it is not necessary to test for the upwind points. The second order algorithm first apply the RT first order algorithm. If the first order phase need to be updated, it means that the wind indeed follows the current sweeping direction and we can then apply the second order correction. It is not a good idea to try to accelerate by minimizing again between the first order and second order phase as this latter can indeed ) alternate directions and their quadrant INRIA be slightly larger. So we just use the second order correction systematically when moving along the wind. There are several subtleties linked to the use of the staggered grid (axes tilted number of points ...). While implementing I significantly raised the number of bugs by insisting in having the same directionnal sweeping direction subroutine and rotating the phase and index of refraction unknown .... We start with the standard homogeneous point sources cases. Phase is prescribed in neighborhood of the sources (several grid points). We use large values at the boundary to emulate outgoing boundary conditions (on a band 2 grid points wide of course ...).

INRIA

Figure 9 and 10 show the contour line for the converged solution for 1 and 2 point sources for :

1. black : first order RT h = 2 200 (10000 grid points). Can be considered as exact. 2. blue : first order RT h = 2 20 . 3. red : second order RT h = 2 20 . As expected second order is much more accuracte and corrects well the anisotropy of the first order scheme along the vertical and horizontal. (linear approximation + curvature), remember we are staggered ... Also the remark that first order performs well along the axis is confirmed. We continue with an heterogeneous test case : (x, y) ∈ [-1, 1] 2 ) n(x, y) = 2.8 if x 2 + y 2 ≥ 0.9 n(x, y) = 2.8 + 20 exp -20 (x 2 +y 2 ) else. It is similar to the heterogenous case of section 5 but the heterogenity is much stronger and rays turn significantly as shown in the contour line plot of figure 14. Initial (boundary) condition is a plane wave (φ = 0) on y = -1. We use the exact boundary conditions φ = 2.8 |y + 1| on x = ±1 and outgoing boundary condition on y = 1. On figures 15 and 16 We plot the phase at the final depth y = 1 for the first order and second order scheme respectively. We use different discretizations N (h = 2/N ).

The first order is clearly first order. The second order scheme converges much faster (it does extermely well when fronts are curved and smooth) except at the kink where it cannot be second order where it still converges without any spurious oscillations thanks to the upwind property of our scheme.

as long as ∂Y ∂y0 does not vanish, Y (., y 0 ) is the ray shot from Γ(y 0 ) and J = det( ∂Y ∂(t,y0) ) is the volume variation of the ray field around this ray. It is therefore necessary to compute J to solve (43). ∂Y ∂t = ∇φ(Y ) is already a computed quantity, we also need ∂Y ∂y0 . A simple Eulerian strategy strategy is to define a "label" variable ξ 0 such that ξ 0 (Y (., y 0 )) = y 0 . Then, a 0 satisfies the even simplest transport equation

∇ξ 0 • ∇φ = 0 ( 45
)
and ∇ξ 0 (Y (., y 0 )) = ( ∂Y ∂y0 ) -1 .

As outlined by the above discussion, second order accuracy is pretty important both for φ as we use the gradient as advecting vector and also because the quantity we are generally interested in is the gradient of ξ 0 . In this section we explain how to adapt the second order scheme of section 5 to solve (45). Extension to a non homogeneous right hand side is straightforward and arises in travel time based tomography [START_REF] Symes | Kirchhoff simulation migration and inversion using finite-difference travel-times and amplitudes[END_REF].

The main idea is still the projection on the coordinate system made of the Lagrangian direction ∇φ and its orthogonal ∇φ ⊥ . We will again use the Lagrangian structure of the Eikonal equation but also the following two fundamental properties of the transport equation (45). First it can be rewritten with as new unknown scalar function α as

∇ξ 0 = α ∇φ ⊥ (46)
Then assuming second order derivative exists (which is necessary to derive a second order scheme anyway) we can take the gradient of (45) which yields

D 2 ξ 0 • ∇φ = -D 2 φ • ∇ξ 0 = -α D 2 φ • ∇φ ⊥ (47)
RR n°0123456789

Figure 2 :

 2 Figure 2: Homogeneous index of refraction : EBRT : red, Exact ray : black -

Figure 3 :

 3 Figure 3: Heterogeneous index of refraction : EBRT : red, Exact ray : black -

Figure 4 :

 4 Figure 4: Authorized angles for the second order scheme

Figure 5 :

 5 Figure 5: Staggered grid and upwind quadrant (shaded)

Figure 6 :

 6 Figure 6: Convex lense

Figure 7 Figure 8 :

 78 Figure 7: Concave lense

Figure 9 :

 9 Figure 9: Homogeneous space : 1 point source

Figure 10 :

 10 Figure 10: Homogeneous space : 2 point sources

Figure 11 :JDBFigure 13 :Figure 14 :

 111314 Figure 11: surface plot of the error. N=20

Figure 15 :

 15 Figure 15: Heterogenous test case : first order scheme -N=10 (green), 20 (blue), 40 (red) -Black : exact solution (first order N=500)

Figure 16 :

 16 Figure 16: Heterogenous test case : second order scheme -N=10 (green), 20 (blue), 40 (red) -Black : exact solution (first order N=500)

  1 N ). Errors are given in thousandth unit (x10 -3 )

	Homogeneous index :
	N	20	40	80	160
	Phi	8.342164822988 4.183548796327 2.096150358148 1.049334449310
	EBRT 7.634700799533 3.671886961461 1.746799748012 0.910111940347
	Heteogeneous index :
	N	20	40	80	160
	Phi	10.141212310403 5.280400011813 2.940170090754 1.485838444092
	EBRT 20.805149751783 8.463070495399 1.555524315021 0.914009038344

  -4 )

	N	10	20	40	80
	Phi	10.28491160694 2.383101764422690 0.5796556601644909 0.1739618008045341
	Phi x 60.08579687216	10.28491160694	2.957285321135705	0.4339775249684319

RR n°0123456789

Acknowledgments : Many thanks to Jianlian Qian, Hong-Kai Zhao, Songtin Luo and Bill Symes for fruitful and helpful discussions. Part of this work was done while visiting Rice university and UC Irvine.

JDB

Using again the decomposition

A taylor second order Taylor expansion gives (using (46) (47)

+O(δX 3 ) (48) In the above equation, the phase gradient ∇φ = (a, b) ( ∇φ ⊥ = (-b, a)) is already given as well as C =< ∇φ ⊥ , D 2 φ • ∇φ ⊥ > the curvature term. They are in principle second order accurate computed quantity of the second order scheme. So the two unknown of our approximation are α and D =< ∇φ ⊥ , D 2 a 0 • ∇φ ⊥ >. Let us write (48) on the same four point stencil ({0, 0}, {-h, 0}, {0, -h}, {-h, -h}) , we revert to the simpler ξ notation instead of ξ 0

(49) Following the upwind relaxation strategy used for the Eikonal equation and after truncating the third order error, the above system must be solved for (ξ 0,0 , α, D) with (ξ 0,-h , ξ -h,0 , ξ -h,-h ) given along with all the "phase" quantities (a, b, C). Looks doable to me ...


Sophie Marat
email: sophie.marat@lis.inpg.frwww.gipsa-lab.inpg.fr

Mickäel Guironnet

Denis Pellerin

VIDEO SUMMARIZATION USING A VISUAL ATTENTION MODEL

This paper presents a method of video summarization based on a visual attention model. The visual attention model is a bottom-up one composed of two parallel ways. A static way, biologically inspired, which highlights salient objects. A dynamic way which gives information about moving objects. A three steps summary method is then presented. The first step is the choice between the two kinds (static and dynamic) of saliency maps given by the attention model. The second step is the selection of keyframes. An "attention variation curve" which highlights changes on frames content during the video is introduced. Keyframes are selected on this variation attention curve. To evaluate the summary a reference summary is built and a comparison method is proposed. The results provide a quantitative analysis and show the efficiency of the video summarization method.

INTRODUCTION

With the development of numerical technologies and Internet the volume of videos has increased spectacularly. Therefore it becomes more and more difficult to retrieve interesting videos or excerpts of video and applications such as video summarization, classification or video browsing are now required. Video summarization goal is to make navigation inside video bases and extraction of important events easier. Two kinds of video summary exist: static ones or dynamic ones, called video skimming. Static video summaries reduce video to few representative frames called keyframes. They are presented like a storyboard, while video skimming is presented like movie previews. In this paper, we consider static video summary which need less storage space.

Most of the video summaries are based on low level features (colour, texture,...) [START_REF] Li | An overview of video abstraction tecnhiques[END_REF] [START_REF] Zhang | An integrated system for content-based video retrieval and browsing I[END_REF] which are not representative of the video semantic content. We describe a method of automatic static video summarization using a higher level feature given a visual attention model. This attention model generates saliency maps, which point out the areas of frames containing more information and attractive for human gaze.

These saliency maps are used to detect changes on frames of the video which make it possible to select keyframes.

The rest of the paper is organised as follows. The attention model used for summarization is presented in section 2. The summarization method is exposed in section 3. In section 4, a summary evaluation method is described. The results are discussed in section 5.

ATTENTION MODEL

The most famous visual attention model is the one proposed by Itti and Koch [START_REF] Itti | A model of saliency-based visual attention for rapid scene analysis[END_REF]. It is a bottom-up (stimulus driven) model based on features like colour, intensity and orientation; it returns a saliency map per frame. A saliency map is a grey scale frame where bright areas that correspond to the regions which attract human gaze. This static model which considers frames one by one, has been improved recently by integrating motion [START_REF] Itti | Realistic avatar eye and head animation using a neurobiological model of visual attention[END_REF].

Ma et al [START_REF] Ma | A generic framework of user attention model and its application in video summarization[END_REF]. have also proposed an user attention model. This model uses a lot of features like static saliency, motion saliency, camera motion, face recognition, aural saliency, keywords, keytopics,... An attention modelling is built on visual attention, aural attention and linguistic attention. This user attention modelling is used for video skimming and a static video can be deduced. This attention model is too complex for static summary only. All the aural features are completely useless for static summary.

The attention model we use is another bottom-up model [START_REF] Guironnet | Static and dynamic feature-based visual attention model: comparison with human judgement[END_REF]. It takes motion into account and is more simpler than the one proposed by Ma. It is built on two parallel ways (figure 1). A static way used to extract textured and contrasted regions of frame. This biologically inspired way gathers a retinal filter, a bank of Gabor filters and interactions between filters answers. A dynamic way is used to detect moving objects by an estimation and a compensation of camera motion. This model gives, after temporal filtering and normalization, a saliency map for each way.

VIDEO SUMMARIZATION METHOD

A video summary from attention model has been proposed by Ma et al. [START_REF] Ma | A generic framework of user attention model and its application in video summarization[END_REF]. It consists in converting the succession of saliency maps into an "attention curve" and selecting keyframes at the maxima of this attention curve. The attention curve is obtained by replacing each saliency maps by the value of the average of its grey levels. The maxima of this curve correspond to the more salient and then contrasted and textured frames. The drawback of this method is that it keeps temporally nearby frames that are likely to have too similar content.

Method principle

In this paper we propose a video summarization method based on the attention model presented in figure 1. It is composed of three steps: the choice of the kind of saliency maps to use (static or dynamic), then the keyframes selection, and finally an elimination of redundant frames.

Choice of the kind of saliency map

This step corresponds to the choice of the more relevant kind of saliency map for each shot of the video to summarize. Indeed, static saliency maps highlight the textured and contrasted objects and are then adapted to describe the shots where motions are weak. At the opposite, dynamic maps highlight moving objects and are adapted for important motions of camera or objects. An attention curve (as defined by Ma et al.) is computed for each kind of saliency map. The attention curve with the higher standard deviation on amplitude corresponds to the kind of map able to give more information and so this is this kind of map that is chosen. Experiments show that this criteria of choice gives the best final results.

Keyframes selection

For the second step, the purpose is to built a summary with a reduce number of frame, but sufficiently different to represent the whole video. For that we propose to generate an "attention variation curve" on which keyframes are selected. To highlight changes during a video, a difference k D of sa- liency map is computed:

k k k i D = M M - -
where k M is the saliency map of the frame k and i the parameter which defines the difference between current map and the i th previous one. Experiments showed that i=10 allows to have frames different enough for seeing changes and to take into account a rapid variation of saliency maps (to lighten when there are changes and to darken when there are no changes).

The attention variation curve is then obtained by computing the average of grey level on each difference map k D .

The selection of keyframes is done by detecting on the attention variation curve (example figure 2), the significant increase with regard to previous neighbours. An adaptive threshold is used for that purpose. The threshold is given by the expression µ +3σ where µ is the average on a sliding window and σ is the standard deviation on the same window. This threshold takes the past into account, and permits to judge if the current frame is different enough from the previous ones. Tests have been made on shots with window's length of 10, 25, 50, 75 and 100 frames. Best results are obtained with a window of 50 frames (i.e. 2 seconds of video). Thus the threshold considers not too much or too few information about the past. A frame is selected each time the variation attention curve crosses the threshold.

Elimination of redundant frames

A summary has been obtained for each shot. The summary was made trying not to forget events in the video and not taking neighbours frames, but frames with similar content can still be selected. A post-processing to eliminate redundant frames is proposed (figure 4). This process considers the static saliency map of each keyframe. A comparison between the saliency maps of the current frame and the previously selected one is done by computing their absolute difference. Then the average of the grey level of this absolute difference is calculated. If this value is below a threshold, the frames are too similar and so only the frames with the highest value on the attention curve are kept. Tests permit to choose the threshold equal to 0.11 and to be sure of its robustness. An example of summary produced by this method is given figure 3. This method is used on long shot but an adaptation has been made for short shots.

Particular case of short shots

The particular case of short shots is considered here. A shot is an unstopped portion of video with or without continuous camera motion. We define a short shot as a shot with 100 frames or less, or 4 seconds of video at a video rate of 25 frames per second. All the frames of a short shot are close temporally. With the continuous variation of frames in video, all the frames of a 4 seconds shot have close content. Only one keyframe is sufficient to summarize that kind of shot. This keyframe is selected by taking the maximum of the attention curve using static saliency map. Static saliency maps are used because in 4 seconds or less the content of the shot do not change a lot and so the static saliency maps which highlight texture and contrast are more efficient than the dynamic ones which highlight motion.

SUMMARY EVALUATION METHOD

Different evaluation method exists. A possibility is to ask subject to choose between two summaries [START_REF] Corchs | Video summarization using a neurodynamical model of visual attention[END_REF]. In this paper the automatic summary is evaluated by comparing it to a "ground-truth" summary called reference summary.

Reference summary

To evaluate the summary previously generated, we compare it to a reference summary. To obtain this summary subjects are asked to watch the shots of a video and to make their ideal summary. They are told to make a summary, with from 1 to 3 keyframes for each shot. For a shot the number N of keyframes for the reference summary is given by the median of the number of frames selected by all the subjects. According to the number of frame chosen by a subject for a given shot a weight is given to each one of these frames. The more frames a subject selects the less their weight are. For each frame, the weight are summed for all the subjects. The reference summary keyframes are selected by taking the N frames with the highest weight.

Comparison method

The comparison between the automatic summary, called candidate summary, and the reference summary is carried out in 4 steps (figure 5).

Figure 1 -

 1 Figure 1 -Principle of attention model.

Figure 2 -

 2 Figure 2 -Application of the redundant frame elimination process to two relevant cases.

Figure 5 -

 5 Figure 5 -three first step of the comparison method

Figure 3 -Figure 4 -

 34 Figure 3 -Saliency maps example for the educative program video. a) frames from the video, b) static saliency maps, c) saliency map difference, d) variation attention curve

Table 1 :

 1 Results for 3 summarization methods on 3 videos (methods: n°1 random summary, n°2 summary selecting one frame at the middle of each shot, n°3 summary with the visual attention model ).

	summary	educative program			TV news			series
		R	P	F1	R	P	F1	R	P	F1
	n°1	62 (15/24) 40 (15/37) 49.1 83 (46/55) 50 (46/91) 63.0 80 (24/30) 40 (24/59) 53.9
	n°2	50 (12/24) 60 (12/20) 54.5 63 (35/55) 83 (35/42) 72.1 73 (22/30) 78 (22/28) 75.8
	n°3	62 (15/24) 51 (15/29) 56.6 78 (43/55) 70 (43/61) 74.1 80 (24/30) 75 (24/32) 77.4

EXPERIMENTAL RESULTS

The summarization method has been tested on 3 videos: an educative program, TV news and a series. All these videos represent 90 shots from 19 to 1468 frames.

The summary evaluation is done using the comparison method describe in section 4 and the criteria of recall R, precision P and their harmonic average F1:

where f N is the number of frames of the candidate sum- mary which correspond with frames of the reference summary, r N the number of frames of the reference summary, and c

N the number of frames of the candidate summary.

Two other methods of summarization (a random summary and a method which selects a keyframe in the middle of each shot) are also compared to the reference summary. Results are given in table 1. For example if we consider the educative program and the summary using attention model (n°3). The recall (respectively the precision) shows that among the 24 (resp. 29) frames of the reference summary the automatic summary retrieves 15 frames.

The summary with visual attention model gives the best results. But the difference with the summary taking a keyframe in the middle of each shot is not high. This is explained by the fact that most of the shots are short shots (less than 100 frames). As we have already seen, the frames of short shots are temporally close and because of the slowness of the variations in video they are semantically close. Thus the two methods give similar results for short shots. The proposed method is more efficient on long shots.

CONCLUSION

In this paper we have presented a summarization method using a visual attention model. This method uses saliency map to compute changes and selects keyframes when they are changes in the video. This summarization method has been tested on three video of different length and content. For that purpose a reference summary and a comparison method were made-up. Results are quite good and improve with the shot length's.

The proposed summarization method has been carried out shot by shot. It could be extended to global summary of the video for reducing the redundancy and being more compact.


Stef Graillat
email: stef.graillat@lip6.fr

Accurate Floating Point Product and Exponentiation

Keywords: accurate product, exponentiation, finite precision, floating point arithmetic, faithful rounding, error-free transformations AMS Subject Classifications: 65-04, 65G20, 65G50

Several different techniques and softwares intend to improve the accuracy of results computed in a fixed finite precision. Here we focus on a method to improve the accuracy of the product of floating point numbers. We show that the computed result is as accurate as if computed in twice the working precision. The algorithm is simple since it only requires addition, subtraction and multiplication of floating point numbers in the same working precision as the given data. Such an algorithm can be useful for example to compute the determinant of a triangular matrix and to evaluate a polynomial when represented by the root product form. It can also be used to compute the power of a floating point number.

Introduction

In this paper, we present fast and accurate algorithms to compute the product of floating point numbers. Our aim is to increase the accuracy at a fixed precision. We show that the results have the same error estimates as if computed in twice the working precision and then rounded to working precision. Then we address the problem on how to compute a faithfully rounded result, that is to say one of the two adjacent floating point numbers of the exact result.

This paper was motived by papers [START_REF] Ogita | Accurate sum and dot product[END_REF][START_REF] Rump | Accurate floating-point summation[END_REF][START_REF] Graillat | Compensated Horner scheme. Research Report 04[END_REF][START_REF] Langlois | How to ensure a faithful polynomial evaluation with the compensated horner algorithm[END_REF] and [START_REF] Kornerup | Computing integer powers in floating-point arithmetic[END_REF] where similar approaches are used to compute summation, dot product, polynomial evaluation and power.

The applications of our algorithms are multiple. One of the examples frequently used in Sterbenz's book [START_REF] Sterbenz | Floating-point computation[END_REF] is the computation of the product of some floating point numbers. Our algorithms can be used to compute the determinant of a triangle matrix

T =      t 11 t 12 • • • t 1n t 22 t 2n . . . . . . t nn      .
Indeed, the determinant of T is

det(T ) = n i=1 t ii .
Another application is for evaluating a polynomial when represented by the root product form p(x) = a n n i=1 (x -x i ). We can also apply our algorithms to compute the power of a floating point number.

The rest of the paper is organized as follows. In Section 2, we recall notations and auxiliary results that will be needed in the sequel. We present the floating point arithmetic and the so-called error-free transformations. In Section 3, we present a classic algorithm to compute the product of floating point numbers. We give an error estimate as well as a validated error bound. We also present a new compensated algorithm together with an error estimate and a validated error bound. We show that under mild assumptions, our algorithm gives a faithfully rounded result. In Section 4, we apply our algorithm to compute the power of a floating point number. We propose two different algorithms: one with our compensated algorithm, the other one with the use of a double-double library.

2 Notation and auxiliary results

Floating point arithmetic

Throughout the paper, we assume to work with a floating point arithmetic adhering to IEEE 754 floating point standard in rounding to nearest [START_REF]IEEE Standard for Binary Floating-Point Arithmetic, ANSI/IEEE Standard[END_REF]. We assume that no overflow nor underflow occurs. The set of floating point numbers is denoted by F, the relative rounding error by eps. For IEEE 754 double precision, we have eps = 2 -53 and for single precision eps = 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. Floating point operations in IEEE 754 satisfy [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF] fl(a

• b) = (a • b)(1 + ε 1 ) = (a • b)/(1 + ε 2 ) for • = {+, -, •, /} and |ε ν | ≤ eps. This implies that |a • b -fl(a • b)| ≤ eps|a • b| and |a • b -fl(a • b)| ≤ eps| fl(a • b)| for • = {+, -, •, /} (2.1)

Error-free transformations

One can notice that a•b ∈ R and fl(a•b) ∈ F but in general we do not have a•b ∈ F. It is known that for the basic operations +, -, •, the approximation error of a floating point operation is still a floating point number (see for example [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]):

x = fl(a ± b) ⇒ a ± b = x + y with y ∈ F, x = fl(a • b) ⇒ a • b = x + y with y ∈ F. (2.2)
These are error-free transformations of the pair (a, b) into the pair (x, y). Fortunately, the quantities x and y in (2.2) can be computed exactly in floating point arithmetic. For the algorithms, we use Matlab-like notations. For addition, we can use the following algorithm by Knuth [START_REF] Knuth | The Art of Computer Programming[END_REF]Thm B. p.236]. Another algorithm to compute an error-free transformation is the following algorithm from Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]. The drawback of this algorithm is that we have x + y = a + b provided that |a| ≥ |b|. Generally, on modern computers, a comparison followed by a branching and 3 operations costs more than 6 operations. As a consequence, TwoSum is generally more efficient than FastTwoSum. For the error-free transformation of a product, we first need to split the input argument into two parts. Let p be given by eps = 2 -p and define s = ⌈p/2⌉. For example, if the working precision is IEEE 754 double precision, then p = 53 and s = 27. The following algorithm by Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] splits a floating point number a ∈ F into two parts x and y such that a = x + y and x and y nonoverlapping with |y| ≤ |x|. With this function, an algorithm from Veltkamp (see [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]) enables to compute an error-free transformation for the product of two floating point numbers. This algorithm returns two floating point numbers x and y such that

a • b = x + y with x = fl(a • b).
Algorithm 2.4 (Veltkamp [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]). Error-free transformation of the product of two floating point numbers function

[x, y] = TwoProduct(a, b) x = fl(a • b) [a 1 , a 2 ] = Split(a) [b 1 , b 2 ] = Split(b) y = fl(a 2 • b 2 -(((x -a 1 • b 1 ) -a 2 • b 1 ) -a 1 • b 2 ))
The following theorem summarizes the properties of algorithms TwoSum and TwoProduct. (2.4)

The algorithm TwoProduct requires 17 flops.

The TwoProduct algorithm can be re-written in a very simple way if a Fused-Multiplyand-Add (FMA) operator is available on the targeted architecture [START_REF] Nievergelt | Scalar fused multiply-add instructions produce floating-point matrix arithmetic provably accurate to the penultimate digit[END_REF]. Some computers have a Fused-Multiply-and-Add (FMA) operation that enables a floating point multiplication followed by an addition to be performed as a single floating point operation. The Intel IA-64 architecture, implemented in the Intel Itanium processor, has an FMA instruction as well as the IBM RS/6000 and the PowerPC before it. On the Itanium processor, the FMA instruction enables a multiplication and an addition to be performed in the same number of cycles than one multiplication or one addition. As a result, it seems to be advantageous for speed as well as for accuracy.

Theoretically, this means that for a, b, c ∈ F, the result of FMA(a, b, c) is the nearest floating point number of a

• b + c ∈ R. The FMA satisfies FMA(a, b, c) = (a • b + c)(1 + ε 1 ) = (a • b + c)/(1 + ε 2 ) with |ε ν | ≤ eps.
Thanks to the FMA, the TwoProduct algorithm can be re-written as follows which costs only 2 flops.

Algorithm 2.5 (Ogita, Rump and Oishi [START_REF] Ogita | Accurate sum and dot product[END_REF]). Error-free transformation of the product of two floating point numbers using an FMA.

function [x, y] = TwoProductFMA(a, b) x = a • b y = FMA(a, b, -x)
3 Accurate floating point product

In this section, we present a new accurate algorithm to compute the product of floating point numbers. In Subsection 3.1, we recall the classic method and we give a theoretical error bound as well as a validated computable error bound. In Subsection 3.2, we present our new algorithm based on a compensated scheme together with a theoretical error bound. In Subsection 3.3, we give sufficient conditions on the number of floating point numbers so as to get a faithfully rounded result. Finally, in Subsection 3.4, we give a validated computable error bound for our new algorithm.

Classic method

The classic method for evaluating a product of n numbers a = (a 1 , a 2 , . . . , a n )

p = n i=1 a i
is the following algorithm.

Algorithm 3.1. Product evaluation

function res = Prod(a) p 1 = a 1 for i = 2 : n p i = fl(p i-1 • a i ) end res = p n
This algorithm requires n -1 flops. Let us now analyse its accuracy. We will use standard notations and standard results for the following error estimations (see [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF]). The quantities γ n are defined as usual [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF] by

γ n := neps 1 -neps for n ∈ N.
When using γ n , we implicitly assume that neps ≤ 1. A forward error bound is

|a 1 a 2 • • • a n -res| = |a 1 a 2 • • • a n -fl(a 1 a 2 • • • a n )| ≤ γ n-1 |a 1 a 2 • • • a n |. (3.5)
Indeed, by induction,

res = fl(a 1 a 2 • • • a n ) = a 1 a 2 • • • a n (1 + ε 2 )(1 + ε 3 ) • • • (1 + ε n ), (3.6
)

with ε i ≤ eps for i = 2 : n. It follows from Lemma 3.1 of [8, p.63] that (1 + ε 2 )(1 + ε 3 ) • • • (1 + ε n ) = 1 + θ n where |θ n-1 | ≤ γ n-1 .
A convenient device for keeping track of power of 1 + ε term is describe in [8, p.68]. The relative error counter k denotes the product

k = k i=1 (1 + ε i ), |ε i | ≤ eps.
A useful rule for the counter is j k = j + k . Using this notation, Equation (3.6) can be written res = fl(a

1 a 2 • • • a n ) = a 1 a 2 • • • a n n -1 .
It is shown in [START_REF] Ogita | Verified solution of linear systems without directed rounding[END_REF] that for a ∈ F, we have

(1 + eps) n ≤ 1 (1 -eps) n ≤ 1 1 -neps , |a| 1 -neps ≤ fl |a| 1 -(n + 1)eps . (3.7)
From Equation (3.6), it follows that

|a 1 a 2 • • • a n -res| ≤ (1 + eps) n-1 γ n-1 |res|.
If meps ≤ 1 for m ∈ N, fl(meps) = meps and fl(1 -meps) = 1 -meps. Therefore,

γ m ≤ (1 + eps) fl(γ m ). (3.8)
Hence,

|a 1 a 2 • • • a n -res| ≤ (1 + eps) n fl(γ n-1 )|res| ≤ (1 + eps) n+1 fl(γ n-1 |res|),
and so

|a 1 a 2 • • • a n -res| ≤ fl γ n-1 |res| 1 -(n + 2)eps .
The previous inequality gives us a validated error bound that can be computed in pure floating point arithmetic in rounding to nearest.

Compensated method

We present hereafter a compensated scheme to evaluate the product of floating point numbers, i.e. the error of individual multiplication is somehow corrected. The technique used here is based on the paper [START_REF] Ogita | Accurate sum and dot product[END_REF]. Algorithm 3.2. Product evaluation with a compensated scheme function res = CompProd(a)

p 1 = a 1 e 1 = 0 for i = 2 : n [p i , π i ] = TwoProduct(p i-1 , a i ) e i = fl(e i-1 a i + π i ) end res = fl(p n + e n )
This algorithm requires 19n -18 flops if we use TwoProduct. It only requires 4n -3 flops if we use TwoProductFMA instead of TwoProduct (if, of course, an FMA is available). Algorithm 3.3. Product evaluation with a compensated scheme with TwoProductFMA

function res = CompProdFMA1(a) p 1 = a 1 e 1 = 0 for i = 2 : n [p i , π i ] = TwoProductFMA(p i-1 , a i ) e i = fl(e i-1 a i + π i ) end res = fl(p n + e n )
It finally costs 3n -2 flops if we also use e i = FMA(e i-1 , a i , π i ) instead of e i = fl(e i-1 a i + π i ).

Algorithm 3.4. Product evaluation with a compensated scheme with TwoProductFMA and FMA

function res = CompProdFMA2(a) p 1 = a 1 e 1 = 0 for i = 2 : n [p i , π i ] = TwoProductFMA(p i-1 , a i ) e i = FMA(e i-1 , a i , π i ) end res = fl(p n + e n )
We will provide an error analysis only for Algorithm CompProd. The error analysis for CompProdFMA1 is the same as CompProd since they share the sames operations. For the error analysis of CompProdFMA2, little changes must be done to take into account the operation e i = FMA(e i-1 , a i , π i ). This changes nearly nothing so it is straightforward to modify the analysis to deal with it.

For error analysis, we note that

p n = fl(a 1 a 2 • • • a n ) and e n = fl n i=2 π i a i+1 • • • a n .
We also have

p = a 1 a 2 . . . a n = fl(a 1 a 2 . . . a n ) + n i=2 π i a i+1 • • • a n = p n + e, (3.9)
where e = n i=2 π i a i+1 • • • a n . Before proving the main theorem, we will need two technical lemmas. The next lemma makes it possible to obtain a bound on the individual error of the multiplication namely π i in function of the inital data a i . Lemma 3.1. Suppose floating point numbers π i ∈ F, 2 ≤ i ≤ n are computed by the following algorithm

p 1 = a 1 for i = 2 : n [p i , π i ] = TwoProduct(p i-1 , a i ) end Then, |π i | ≤ eps(1 + γ i-1 )|a 1 • • • a i | for i = 2 : n.
Proof. From Equation (2.1), it follows that

|π i | ≤ eps|p i |.
Moreover, p i = fl(a 1 • • • a i ) so that from (3.5),

|p i | ≤ (1 + γ i-1 )|a 1 • • • a i |.
Hence,

|π i | ≤ eps(1 + γ i-1 )|a 1 • • • a i |.
The following lemma enables us to bound the rounding errors during the computation of the error during the full product. Lemma 3.2. Suppose floating point numbers e i ∈ F, 1 ≤ i ≤ n are computed by the following algorithm

e 1 = 0 for i = 2 : n [p i , π i ] = TwoProduct(p i-1 , a i ) e i = fl(e i-1 a i + π i ) end Then, |e n - n i=2 π i a i+1 • • • a n | ≤ γ n-1 γ 2n |a 1 a 2 • • • a n |.
Proof. First, one notices that e n = fl( n i=2 (π i a i+1 • • • a n )). We will use the error counters described above. For n floating point numbers x i , it is easy to see that [8, chap.4] fl(

x 1 + x 2 + • • • + x n ) = x 1 n -1 + x 2 n -1 + x 3 n -2 + • • • + x n 1 .

This implies that

e n = fl( n i=2 (π i a i+1 • • • a n )) = fl(π 2 a 3 • • • a n ) n -2 + fl(π 3 a 4 • • • a n ) n -2 + • • • + fl(π n ) 1 .
Furthermore, we have shown before that fl(a

1 a 2 • • • a n ) = a 1 a 2 • • • a n n -1 . Consequently, e n = π 2 a 3 • • • a n n -2 n -1 + π 3 a 4 • • • a n n -3 n -1 + • • • + π n 1 . A straightforward computation yields |e n - n i=2 π i a i+1 • • • a n | ≤ γ 2n-3 n i=2 |π i a i+1 • • • a n |. From Lemma 3.1, we have |π i | ≤ eps(1 + γ i-1 )|a 1 • • • a i | and hence |e n - n i=2 π i a i+1 • • • a n | ≤ (n -1)eps(1 + γ n-1 )γ 2n-3 |a 1 a 2 • • • a n |.
Since eps(1 + γ n-1 ) = γ n-1 /(n -1) and γ 2n-3 ≤ γ 2n , we obtain the desired result.

One may notice that the computation of e n is similar to the Horner scheme. One could have directly applied a result on the error of the Horner scheme [8, Eq.(5.3),p.95].

We can finally state the main theorem.

Theorem 3.3. Suppose Algorithm 3.2 is applied to floating point number a i ∈ F, 1 ≤ i ≤ n, and set p = n i=1 a i . Then, |res -p| ≤ eps|p| + γ n γ 2n |p|.

Proof. The fact that res = fl(p n + e n ) implies that res

= (1 + ε)(p n + e n ) with ε ≤ eps. So it follows |res -p| = | fl(p n + e n ) -p| = |(1 + ε)(p n + e n -p) + εp| = |(1 + ε)(p n + n i=2 π i a i+1 • • • a n -p) + (1 + ε)(e n - n i=2 π i a i+1 • • • a n ) + εp| = |(1 + ε)(e n - n i=2 π i a i+1 • • • a n ) + εp| by (3.9) ≤ eps|p| + (1 + eps)|e n - n i=2 π i a i+1 • • • a n | ≤ eps|p| + (1 + eps)γ n-1 γ 2n |a 1 a 2 • • • a n |. Since (1 + eps)γ n-1 ≤ γ n , it follows that |res -p| ≤ eps|p| + γ n γ 2n |p|.
It may be interesting to study the condition number of the product evaluation. Ones defines

cond(a) = lim ε→0 sup |(a 1 + ∆a 1 )(a 2 + ∆a 2 ) • • • (a n + ∆a n ) -a 1 a 2 • • • a n | ε|a 1 a 2 • • • a n | : |∆a i | ≤ ε|a i | .
A standard computation yields cond(a) = n.

Faithful rounding

We define the floating point predecessor and successor of a real number r satisfying min{f : f ∈ R} < r < max{f : f ∈ F} by pred(r) := max{f ∈ F : f < r} and succ(r) := min{f ∈ F : r < f }.

Definition 3.1. A floating point number f ∈ F is called a faithful rounding of a real number r ∈ R if pred(f ) < r < succ(f ).
We denote this by f ∈ (r). For r ∈ F, this implies that f = r.

A faithful rounding is then one of the two adjacent floating point numbers of the exact result.

that is γ n γ 2n < 1 -eps 2 + eps eps. Since γ n γ 2n ≤ 2(neps) 2 /(1 -2neps) 2 , a sufficient condition is 2 (neps) 2 (1 -2neps) 2 < 1 -eps 2 + eps eps which is equivalent to neps 1 -2neps < (1 -eps)eps 2(2 + eps)
and then to

n < √ 1 -eps √ 2 √ 2 + eps + 2 (1 -eps)eps eps -1/2 .
We have just shown that if n < αeps -1/2 where α ≈ 1/2 then the result is faithfully rounded. More precisely, in double precision where eps = 2 -53 , if n < 2 25 ≈ 5 • 10 7 , we get a faithfully rounded result.

Validated error bound

We present here how to compute a valid error bound in pure floating point arithmetic in rounding to nearest. It holds that

|res -p| = | fl(p n + e n ) -p| = | fl(p n + e n ) -(p n + e n ) + (p n + e n ) -p| ≤ eps|res| + |p n + e n -p| ≤ eps|res| + |e n -e|. Since |e n -e| ≤ γ n-1 γ 2n |p| and |p| ≤ (1 + eps) n-1 fl(|a 1 a 2 • • • a n |) we obtain |res -p| ≤ eps|res| + γ n-1 γ 2n |p| ≤ eps|res| + γ n-1 γ 2n (1 + eps) n-1 fl(|a 1 a 2 • • • a n |).
Using (3.7) and (3.8), we get

|res -p| ≤ fl(eps|res|) + (1 + eps) n fl(γ n ) fl(γ 2n ) fl(|a 1 a 2 • • • a n |) ≤ fl(eps|res|) + (1 + eps) n+2 fl(γ n γ 2n |a 1 a 2 • • • a n |) ≤ fl(eps|res|) + fl γ n γ 2n |a 1 a 2 • • • a n | 1 -(n + 3)eps ≤ (1 + eps) fl eps|res| + γ n γ 2n |a 1 a 2 • • • a n | 1 -(n + 3)eps ≤ fl eps|res| + γ n γ 2n |a 1 a 2 • • • a n | 1 -(n + 3)eps / (1 -2eps) .
We can summarize this as follows.

Lemma 3.7. Suppose Algorithm 3.2 is applied to floating point numbers a i ∈ F, 1 ≤ i ≤ n and set p = n i=1 a i . Then, the absolute forward error affecting the product is bounded according to

|res -p| ≤ fl eps|res| + γ n γ 2n |a 1 a 2 • • • a n | 1 -(n + 3)eps / (1 -2eps) .

Validated error bound and faithful rounding

In the previous subsection, we have shown that

|e n -e| ≤ fl γ n γ 2n |a 1 a 2 • • • a n | 1 -(n + 3)eps . ( 3
γ n γ 2n |a 1 a 2 • • • a n | 1 -(n + 3)eps < fl(eps|res|)
then we got a faitfully rounded result. This makes it possible to check a posteriori if the result is faithfully rounded.

Exponentiation

In this section, we study two exponentiation algorithms (for computing x n with x ∈ F and n ∈ N). The first one is linear (in O(n)) whereas the second one is logarithmic (in O(log n)).

A linear algorithm

The natural method to compute x n is to apply algorithm CompProd for a i ∈ F with a i = x for 1 ≤ i ≤ n. This leads to the following algorithm. All that has been done concerning CompProd is still valid with CompLinPower. For example, if n < 2 25 then the result is faithfully rounded.

This algorithm is similar to the one of [START_REF] Kornerup | Computing integer powers in floating-point arithmetic[END_REF]. It is also the same as the Compensated Horner scheme [START_REF] Graillat | Compensated Horner scheme. Research Report 04[END_REF] applied to the polynomial p(x) = x n . Results concerning faithful polynomial evaluation can be found in [START_REF] Langlois | How to ensure a faithful polynomial evaluation with the compensated horner algorithm[END_REF].

A double-double library

Compensated methods are a possible way to improve the accuracy. Another possibility is to increase the working precision. For this purpose, some multiprecision libraries have been developed. One can divide those libraries into three categories.

• Arbitrary precision library using a multiple-digit format where a number is expressed as a sequence of digits coupled with a single exponent. Examples of this format are Bailey's MPFUN [START_REF] Bailey | A Fortran 90-based multiprecision system[END_REF], Brent's MP [START_REF] Brent | A Fortran multiple-precision arithmetic package[END_REF] or MPFR [START_REF]MPFR, the Multiprecision Precision Floating Point Reliable library[END_REF].

• Arbitrary precision library using a multiple-component format where a number is expressed as unevaluated sums of ordinary floating point words. Examples of this format are Priest [START_REF] Douglas | On Properties of Floating Point Arithmetics: Numerical Stability and the Cost of Accurate Computations[END_REF] and Shewchuk [START_REF] Richard | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF].

• Extended fixed precision library using the multiple-component format but with a limited number of components. Examples of this format are Bailey's double-double [START_REF] Bailey | A Fortran-90 double-double library[END_REF] (doubledouble numbers are represented as an unevaluated sum of a leading double and a trailing double) and quad-double [START_REF] Hida | Algorithms for quad-double precision floating point arithmetic[END_REF] (quad-double numbers are represented as an unevaluated sum of four IEEE doubles).

In the sequel, we present two algorithms to compute product of two double-double or a double times a double-double. Those algorithms are taken from [START_REF] Quirin | Basic building blocks for a triple-double intermediate format[END_REF].

Algorithm 4.2. Multiplication of two double-double numbers function

[r h , r l ] = prod dd dd(a h , a l , b h , b l ) [t 1 , t 2 ] = TwoProduct(a h , b h ) t 3 = fl(((a h • b l ) + (a l • b h )) + t 2 ) [r h , r l ] = TwoProduct(t
] = prod dd d(a, b h , b l ) [t 1 , t 2 ] = TwoProduct(a, b h ) t 3 = fl((a • b l ) + t 2 ) [r h , r l ] = TwoProduct(t 1 , t 3 )
Theorem 4.1 (Lauter [13, thm. 4.7]). Let be a h + a l and b h + b l the double-double arguments of Algorithm 4.2. Then the returned values r h and r l satisfy

r h + r l = ((a h + a l ) • (b h + b l ))(1 + ε)
where ε is bounded as follows : |ε| ≤ 16eps 2 . Furthermore, we have |r l | ≤ eps|r h |.

Results for Algorithm 4.3 are very similar with a = a h and a l = 0.

A logarithmic algorithm

A logarithmic algorithm was introduced in [START_REF] Kornerup | Computing integer powers in floating-point arithmetic[END_REF] using the classic right-to-left binary exponentiation algorithm and the double-double library. Hereafter, we propose a variant of this algorithm with the left-to-right binary exponentiation algorithm together with the double-double library. Contrary to the right-to-left binary exponentiation algorithm which needs two multiplications of two double-double numbers, the left-to-right binary exponentiation algorithm only needs a multiplication of two double-double numbers and a multiplication of a double number by a double-double numbers. Moreover, the multiplication of a double-double by a double-double is actually a square so that it can be a little bit optimized.

h + l = x n (1 + ε) with (1 -16eps 2 ) n-1 ≤ 1 + ε ≤ (1 + 16eps 2 ) n-1 .
Proof. The proof is very similar to the one of [START_REF] Kornerup | Computing integer powers in floating-point arithmetic[END_REF]Thm 4]. It comes from the fact that by induction one can show that the approximation of x k is of the form

x k (1 + ε k ) with (1 - 16eps 2 ) k-1 ≤ 1 + ε k ≤ (1 + 16eps 2 ) k-1 .
Let us denote ϕ = eps(1eps) and

γ n = 16neps 2 1 -16neps 2 .
It follows that 1 + ε ≤ (1 + 16eps 2 ) n-1 ≤ 1 + γ n [8, p.63]. Since h + l = x n (1 + ε) then it holds x n = h + l -εx n . Furthermore, it also holds res = fl(h + l). From Lemma 3.5, if we prove that 2ε|x n | < eps|res| then res is a faithful rounding of x n . From the definition of res, it follows that (1eps)|h + l| ≤ |res| and so (1eps)(1 + ε)|x n | ≤ |res|. A sufficient condition to ensure a faithful rounding is then 2ε|x n | < eps(1eps)(1 + ε)|x n | and so 2ε < eps(1eps)(1 + ε), which is equivalent to ε < ϕ 2 -ϕ .

Since ε ≤ γ n , a sufficient condition is

γ n < ϕ 2 -ϕ ,
which is equivalent to

n < ϕ 16eps 2 =
(1eps)eps -1 16 .

For example, in double precision where eps = 2 -53 , if n < 2 49 ≈ 5 • 10 14 , then we get a faithfully rounded result.

Conclusion

In this paper, we provided an accurate algorithm for computing product of floating point numbers. We gave some sufficient conditions to obtain a faithfully rounded result as well as validated error bounds. We applied this algorithm to compute exponentiation of floating point numbers. We improved this algorithm by using a double-double library.

Algorithm 2 . 1 (

 21 Knuth [10]). Error-free transformation of the sum of two floating point numbers function [x, y] = TwoSum(a, b) x = fl(a + b) z = fl(x -a) y = fl((a -(x -z)) + (b -z))

Algorithm 2 . 2 (

 22 Dekker [5]). Error-free transformation of the sum of two floating point numbers. function [x, y] = FastTwoSum(a, b) x = fl(a + b) y = fl((a -x) + b)

Algorithm 2 . 3 (

 23 Dekker [5]). Error-free split of a floating point number into two parts function [x, y] = Split(a, b) factor = fl(2 s + 1) c = fl(factor • a) x = fl(c -(c -a)) y = fl(a -x)

Theorem 2 . 1 (. 3 )

 213 Ogita, Rump and Oishi [15]). Let a, b ∈ F and let x, y ∈ F such that [x, y] = TwoSum(a, b) (Algorithm 2.1). Then, a + b = x + y, x = fl(a + b), |y| ≤ eps|x|, |y| ≤ eps|a + b|. (2The algorithm TwoSum requires 6 flops. Let a, b ∈ F and let x, y ∈ F such that [x, y] = TwoProduct(a, b) (Algorithm 2.4). Then, a • b = x + y, x = fl(a • b), |y| ≤ eps|x|, |y| ≤ eps|a • b|.

Corollary 3 . 4 .

 34 Suppose Algorithm 3.2 is applied to floating point number a i ∈ F, 1 ≤ i ≤ n, and set p = n i=1 a i . Then, |res -p| |p| ≤ eps + γ n γ 2n n cond(a).

Lemma 3 . 5 (

 35 Rump, Ogita and Oishi [18, lem. 2.5]). Let r, δ ∈ R and r := fl(r). Suppose that 2|δ| < eps| r|. Then r ∈ (r + δ), that means r is a faithful rounding of r + δ.Let res be the result of CompProd. Then we have p = p n + e and res = fl(p n + e n ) with e = n i=2 π i a i+1 • • • a n . It follows that p = (p n + e n ) + (e -e n). This leads to the following lemma which gives a criterion to ensure that the result of CompProd is faithfully rounded. Lemma 3.6. With the previous notations, if 2|e-e n | < eps|res| then res is a faithful rounding of p. Since we have |e -e n | ≤ γ n γ 2n |p| and (1eps)|p| -γ n γ 2n |p| ≤ |res|, a sufficient condition to ensure a faithful rounding is 2γ n γ 2n |p| < eps((1eps)|p| -γ n γ 2n |p|)

Algorithm 4 . 1 .

 41 Power evaluation with a compensated schemefunction res = CompLinPower(x, n) p 1 = x e 1 = 0 for i = 2 : n [p i , π i ] = TwoProduct(p i-1 , x) e i = fl(e i-1 x + π i ) end res = fl(p n + e n )

Algorithm 4 . 4 .

 44 Power evaluation with a compensated schemefunction res = CompLogPower(x, n) % n = (n t n t-1 • • • n 1 n 0 ) 2 [h, l] = [1, 0] for i = t : -1 : 0 [h, l] = prod dd dd(h, l, h, l) if n i = 1 [h, l] = prod dd d(x, h, l) end end res = fl(h + l)Theorem 4.2. The two values h and l returned by Algorithm 4.4 satisfy

  .10) Lemma 3.6 tells us that if 2|e -e n | < eps|res| then res is a faithful rounding of p (where res is the result of CompProd).

	As a consequence, if
	fl 2