H Escamilla Núñez

F Mora Camino

Generation of Optimal Smooth Trajectories using Bezier Curves for Transport Aircraft

Keywords: Flight trajectory, trajectory generation, Bezier curves, optimization, 4D guidance

With the increase in air traffic, new concepts have been introduced to organize and better manage air traffic flows (free flight, air-corridors, airstreams) with the view of increasing traffic safety and airspace capacity by solving air traffic conflicts. Hence the need of ad-hoc descriptions and parametrization of more complex and flexible transport aircraft trajectories with new characteristics, allowing high traffic densities as well as limiting environmental impact (noise). In the first part of this paper the motivation and trajectory generation problem for transportation aircraft is introduced, and a state of the art with respect to transportation aircraft trajectory generation techniques is proposed. Then a smooth path generation based on Bezier curves is proposed.

I. INTRODUCTION

Nowadays, commercial airliners aircraft even with operating Airborne Collision Avoidance Systems (ACAS) or Traffic Alert and Collision Avoidance Systems (TCAS), are exposed to near mid air collisions, or in worst case scenarios, to events like the one occurred on july 1 over Überlingen, Germany in 2002, where two transport-category aircraft collided after TCAS instructed one pilot to climb, but the pilot descent in compliance to air traffic controller instructions. Moreover, even if only 1 in 10 Traffic Advisory (TAs) provided by TCAS result in a Resolution Advisory (RAs), situations like:

• January 2001, a Boeing 747 followed the instruction of the ATC to descend instead of his TCAS RA to climb, intersecting course with a DC-10. The collision was avoided when the 747 was put into a steep descent after visual contact with the other aircraft. However, about 100 crew and passengers aboard the 747 sustained injuries due the emergency maneuver. • Switzerland, June 2011, Raytheon 390 crew followed an ATC descent clearance during their TCAS climb RA, creating a conflict with an Airbus 319. Both aircraft passed in very close proximity (0.6 nm horizontally and 50 feet vertically) without either sighting the other. • Switzerland, May 2012, an Airbus A320 departing Zurich in a climbing turn received a TCAS RA climb caused by an AW 139 only equipped with TCAS 1 also departing Zurich. The conflict in Class 'C' airspace was attributed to inappropriate clearance issue by the TWR controller and their inappropriate separation monitoring. • May 2013, an A319 in Swiss Class 'C' airspace received a TCAS 'Level Off' RA against a 737 above after being inadvertently given an incorrect climb clearance by ATC. and other1 , are clear examples of human error and lack of homogeneity in TCAS. In spite of the fact that up-to-date TCAS (III and IV) versions were capable to give Vertical and/or Horizontal direction as RAs, both were abandoned due to new arising issues. As both SESAR and NextGen projects plan to implement new operational concepts which will reduce the spacing between aircraft. TCAS in its current form is not compatible with such concepts and would alert too frequently to be useful. Consequently, in order to allow aircraft to fly closer, quick and efficient 4D trajectory generation is crucial. In this manner, as position and time of aircraft during fly is properly handled, air traffic conflicts can be avoided or addressed if present, taking into account load factor limits to look after passengers well-being.

The generation of flyable and efficient trajectories has been considered by several authors [START_REF] Bakolas | Initial guess generation for aircraft landing trajectory optimization[END_REF], [START_REF] Anderson | Real-time dynamic trajectory smoothing for unmanned air vehicles[END_REF], [START_REF] Judd | Spline based path planning for unmanned air vehicles[END_REF], [START_REF] Yang | Generation of dynamically feasible and collision free trajectory by applying six order bezier curve and local optimal reshaping[END_REF], [START_REF] Yang | An analytical continuous-curvature path-smoothing algorithm[END_REF], [START_REF] Delahaye | Mathematical Models for Aircraft Trajectory Design: A survey[END_REF]. From a general view, reference trajectories are essential for flight plans which meet a large set of overfly or profile constraints, which vary in general from flight to flight. On the other hand, the generation of reference trajectories at short term should allow the implementation of procedures in the case of potential conflicts. Current path generation for transport aircraft is based on a sequence of objective points in a 2D or 3D space given to connect two geographical locations [START_REF] Walter | Digital Avionics Handbook[END_REF]. Then, using diagrams like Voronoi, or Delaunay triangulation, or any other method, a piecewise path is constructed using straight lines denominated legs. This path may be offered by a high-level path planner from techniques such as Dijkstra's, A*, probabilistic roadmaps, genetic algorithms [START_REF] Sahingoz | Generation of bezier curve-based flyable trajectories for multi-uav systems with parallel genetic algorithm[END_REF], or Rapidly exploring random tree star (RRT*) [START_REF] Webb | Kinodynamic rrt*: Asymptotically optimal motion planning for robots with linear dynamics[END_REF], just to mention a few. However, as path differential requirements have to be satisfied, the initial trajectory needs to be reshaped in order to provide a flyable trajectory for transport aircraft. The problem to design a flyable path over or close the objective points while satisfying constraints such as maximum curvature and/or G 2 continuity (curvature continuity), has been addressed using different approaches. As a pioneer, Dubins assured that the shortest path between two points involves circles and straight line path segments [START_REF] Dubins | On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents[END_REF]. However, curvature continuity at the joints of lines and circular arcs is not satisfied. Techniques to solve this track transition problem are addressed using Clothoids [START_REF] Scheuer | Continuous-curvature path planning for car-like vehicles[END_REF], but as they don't have closed-form expressions, the computation complexity is increased. Circular arcs as transitions tracks have been also proposed [START_REF] Anderson | Real-time dynamic trajectory smoothing for unmanned air vehicles[END_REF], where a 2D real-time trajectory is generated satisfying curvature and velocity constraints. Also, the deviation between the generated circular path and the associated control points is minimized using a parameter κ ∈ [0, 1]. If κ = 0, the distance is zero, so the plane will fly exactly over the waypoint, and if κ = 1, minimum-time transitions between control points are achieved. However, the election of this parameter κ becomes an issue when a range of distance wants to be elected as a permitted deviation from the generated path. Furthermore, Dubins like approaches generate paths limited to straight lines and arcs of circles without parametrization, thus, the method shows disadvantages being able to generate several arcs without curvature continuity or with high computation complexity. Another tool to generate smooth flyable paths are splines, defined by series of low order polynomials. In [START_REF] Judd | Spline based path planning for unmanned air vehicles[END_REF], after finding the best path from a UAV position location to a target location using a Voronoi diagram and a Dijkstra's algorithm, the 2D path is smoothed using cubic splines. It is worth to say that the optimal locations of the middle knots of a spline are crucial for the shape of the segments. However, the optimal computation of these points is very time-consuming unless a set of cases are defined apriori. As an attractive approach, Bezier curves are able to generate continuous-curvature paths, having the advantage of passing through initial and final points while the whole curve always lies within the convex hull that is constructed by the control points. Also, the initial and final piecewise straight lines of the control polygon are always tangent to the Bezier curve at the starting and ending control points. An example is given in [START_REF] Lee | Planning and control for collision-free cooperative aerial transportation[END_REF], where the authors present the generation of a 3D path obtained from a combination of Rapidly exploring random tree star (RRT*) using a C 2 class Bezier curve, and Dynamic Movement Primitives (DMP's), that allow cooperative aerial manipulators to avoid known and unknown obstacles. As the final Bezier curve is a barycentric combination of the polygon vertices formed by the control points, one may think that a trajectory cannot be accurately generated if the control points positions are zig-zags like, however, this aspect can be handled by stitching several Bezier curves to form a bigger path while assuring G 2 continuity at the joints. In [START_REF] Yang | Generation of dynamically feasible and collision free trajectory by applying six order bezier curve and local optimal reshaping[END_REF], after a Guiding Attraction based Random Tree (GART) is used to obtain control points in a 2D plane, a kinematic smoother based on sixth order Bezier curves to achieve second derivative continuity (curvature continuity) is proposed. Finally, a local optimal reshaping of the path, minimizing length and curvature cost is performed. In [START_REF] Yang | An analytical continuous-curvature path-smoothing algorithm[END_REF], it is presented an algorithm based on cubic Bezier curves for 3D path smoothing, satisfying G 2 continuity and maximum curvature constraints, where instead of addressing a direct solution for the 3D path smoothing, a 2D path smoothing for consecutive triplets of control points is applied, seeing each triplet as a 2D plane thanks to the Frenet frame; Nevertheless, deviation from the resulting trajectory w.r.t. the control points is not controlled at all. In this paper, assuming velocity or time constraints, a time-parametrized smooth trajectory valid for 4D guidance is developed by stitching several Bezier curves while assuring G 2 continuity at the joints, also, the Euclidian distance between a certain control point and the proposed trajectory is controlled, yielding an optimal 4D path as a reference trajectory. The paper is organized as follows, Section II gives some theory on Bezier curves and conditions for G 2 continuity. Then the path generation is explained in Section III. In Section IV, the optimal reshaping of the path is described. Also, the obtention of the load factor for independent vertical or lateral trajectories is obtained in Section V Finally, conclusions are given in Section VI.

II. BEZIER CURVES

A Bezier curve P (s) of degree n, obtained from n + 1 control points (P 0 , . . . , P n ), is described by

P (s) = n i=0 B n i (s)P i s ∈ [0, 1] (1)
where B n i (s) is the i th Bernstein polynomial of degree n, given by

B n i (s) = n i s i (1 -s) n-i i ∈ {0, 1, . . . , n} (2)
Let a Bezier curve C(s) with m + 1 control points (C 0 , . . . , C m ), and a second Bezier curve D(s) with n + 1 control points (D 0 , . . . , D n ) be joined. According to [START_REF] Barsky | Geometric continuity of parametric curves[END_REF] on C n class and regularity, C n Continuity, and G n Continuity, C 0 continuity is guaranteed if:

C m = D 0 (3)
Then, knowing that the first curve is tangent to the last leg, and the second curve is tangent to the first leg, a smooth transition is assured if ( 3) is satisfied and C m , D 0 , C m-1 , D 1 are on the same line. Furthermore, C 1 continuity is guaranteed if the tangent vector of the first curve at s = 1 is identical to the tangent vector of the second curve at s = 0, meaning that:

C (1) = m(C m -C m-1 ) = n(D 1 -D 0 ) = D (0) (4)
This states that the ratio Cm-C m-1

D 1 -D 0
, involving the length of the last leg of the first curve ( C m -C m-1 ), and the length of the first leg of the second curve ( D 0 -D 1 ), must be n m . Since n and m are fixed numbers, the positions of C m-1 and D 1 can be rearranged to be not only at the same line, but also at the proper distance to assure C 1 continuity. In the same tenor, C 2 and G 2 continuity are guaranteed at the joint if it is verified that

C (1) = m(m -1)(C m -2C m-1 + C m-2 ) = n(n -1)(D 2 -2D 1 + D 0 ) = D (0) (5)
Hence, to assure G 2 continuity, the locations of

C m , C m-1 , C m-2 , D 0 , D 1 , D 2
, where (3)-( 5) are satisfied, need to be proposed.

III. TRAJECTORY GENERATION A. G 2 Continuity Path

In order to handle the curvature while the torsion (what distinguishes a circle from a helix) is zero, the 3D path is decomposed into several 2D planar trajectories laying on the tangent/normal plane of the frenet frame, as some authors have proposed [START_REF] Yang | An analytical continuous-curvature path-smoothing algorithm[END_REF]. Since three control points are required to form a plane, at least n = 3 control points (P 1 , . . . , P n ) are defined, and then divided into n -2 triplets. Also, in order to assure an accurate generation of the path even if the objective points are zig-zags like, several quintic Bezier curves, one for each triplet of control points, are stitched together to form a bigger path while assuring G 2 continuity at the joints. For a triplet of control points (P 1 , P 2 , P 3 ), like the ones depicted in figure 1, G 1 continuity can be achieved by interpolating four points (Q 0 , Q 1 , Q 4 , Q 5 ), and for G 2 continuity, six points are interpolated (Q i ; i ∈ {0, . . . , 5}). After the interpolation of these auxiliary control points is done, in the case of the G 2 continuity path, a quintic Bezier curve is adjusted to them. The points are interpolated as follows. Q 0 and Q 5 are defined to be at the middle point of ( P1 P 2 ) and ( P2 P 3 ) respectively, forcing a past and following Bezier curve formed by the past and next triplet of control points, to finish and start at Q 0 and Q 5 control points, in other words, a next triplet of points conformed by (P 2 , P 3 , P 4 ) will have a Bezier curve starting at Q 5 , so (3) is fulfilled. For Q 1 and Q 4 , they are computed to be also in the same line of ( P1 P 2 ) and ( P2 P 3 ) respectively, but separated by a δ 1 distance from Q 0 in the case of Q 1 , and a distance δ 2 from Q 5 in the case of Q 4 . Note that the number of auxiliary interpolated points to compute the Bezier curves are the same for all the triplet of points, so m and n from (4) are equal, hence, to fulfill G 1 continuity, only remains to guarantee that the distance δ 2 of a certain triplet of control points is equal to the distance δ 1 of the next triplet of control points. Finally, if the positions of Q 2 and Q 3 are computed also in the same line of ( P1 P 2 ) and ( P2 P 3 ) respectively, separated by the same distances δ 1 from Q 1 , and δ 2 from Q 4 , if it is proposed: δ 1 = P1 P 2 4 , and δ 2 = P2 P 3 4 . Equations ( 4), ( 5) are satisfied, assuring G 2 continuity as the curvature is forced to be zero at the joints between Bezier curves. Besides, knowing that the first and last Bezier curves conforming the total path will start and end at the middle point of their correspondent control points, as straight lines have zero curvature, the path can be completed with straight lines without affecting the G 2 continuity (see figure 1).

B. Parametrization of the path

As several Bezier curves conform the total path, and each Bezier curve is parametrized by s ∈ [0, 1], a time re-parametrization needs to be done for the creation of a flyable path for transport aircraft. First, the arc length of each Bezier curve is obtained by integrating from zero to one the norm of the first derivative of the Bezier curves, as well as the initial and final straight lines that conform the total path. Once the arc lengths are obtained, a velocity for the aircraft to follow the path can be chosen, or as an alternative way, the time in which the aircraft is supposed to fly over the joints of the Bezier curves is defined. For the numerical simulation, some control points were chosen, given in Table 1.

These points generate a G 1 and G 2 path depicted in figure 2. Regarding to the joints of the G 2 Bezier curves, and assuming that the path is intended to be followed at a constant 200m/s velocity, figure 2 shows the arc lengths (l i ; i ∈ [START_REF] Bakolas | Initial guess generation for aircraft landing trajectory optimization[END_REF][START_REF] Delahaye | Mathematical Models for Aircraft Trajectory Design: A survey[END_REF]) and their times (t j ; j ∈ [0, 6]) assigned. Note that these values of time, are the times at which the Bezier curves are joined. The values of l i and t j are given in Table 2. However, as Bezier curves are not parametrized by time, an operation is performed such that:

s = t -t i t i+1 -t i i ∈ [0, 1, ..., 5] (6)
So that the multiple Bezier curves parametrized by s ∈ [0, 1], can be used at the proper time intervals [t i , t i+1 ]. Finally, comparing the curvatures of the G 1 and G 2 continuity paths (see figure 3), it is clearly seen how the G 2 path has continuity in the curvature at the joints of the different Bezier curves, and that these joints occur at the assigned time, represented by asterisks in the correspondent axis. IV. OPTIMAL RESHAPING Depending on the application, some aircraft may need to fly directly over the control points or at least near a defined range of distance from the control points. In order to achieve this demand, the Euclidian distance from the G 2 continuity path w.r.t. the nearest control point is controlled. For this special case, an extra auxiliary point (Q 3 in figure 4) is computed such that an optimum path based on sixth order Bezier curves is obtained. The other auxiliary control points are conserved as they are in order to do not affect curvature continuity, but the position of the heptic auxiliary point will be moved gradually until the generated path passes within a distance range defined by the user. The direction in which this point is moved will be denoted by

aux = (Q 2 -Q 1 ) + (Q 4 -Q 5 ) (Q 2 -Q 1 ) + (Q 4 -Q 5 ) (7)
In figure 4, a 100m maximum deviation for a triplet of control points is commanded, meaning that the optimized path will be generated such that it will pass no further than 100m away the control point P 2 . In this case, the optimizing algorithm stopped at a distance of 92.48m. The distance at which the algorithm stops depends on the step size in which Q 3 is moved away in the direction of aux. For a small step size, the path will be closer to the distance defined by the user, in this case 100m, but the computation time will be increased. Now, for a complete trajectory, using the control points of Table 1, a 100m maximum deviation is commanded, the resultant path is shown in figure 5.

The Euclidian distance at which the trajectory is generated from the closest control points is given by Table 3. Assuring that the maximum deviation distance is achieved. As close as the optimal path is to the control points, the curvature becomes bigger, such that a maximum curvature could also be defined and controlled by the deviation distance defined by the user, establishing a tradeoff between deviation of the path from the control points and maximum curvature constraints. Furthermore, as the generated trajectory is time-parametrized, it could be used for 4D guidance problems.

V. LOAD FACTOR OF THE TRAJECTORY

As transport aircraft are designed to flight smooth trajectories, some limitations on the load factor have to be respected for the well being and comfort of passengers. According to [14], on a typical flight, the load factor is limited to +2.5g and -1g for regular maneuvers, or up to 3.8g in maximum takeoff weight. The expressions for the load factor (in the body frame) of an aircraft are given by Under the assumptions of a steady turn, it can be obtained that

n x = -s θ + (rv -qw)/g (8) n y = c θ s φ + (pw -ru)/g (9) n z = c θ c φ + (qu -pv)/g (10)
n z = c θ c φ (11)
Note that the bigger value of n z is obtained when θ = 0 • , so, for a 2.5g limit, the maximum bank angle is 66.5 • . Also, an analysis involving the centrifugal acceleration yields

mV 2 a R = Ls φ (12a) Lc φ = W (12b)
thus, (12b) can be rewritten as

L W = 1 c φ = n z (13)
which is also when θ = 0 • in (11), then, using this relation into (12a), the load factor can be related with the radius of horizontal curvature R as

R = 1 n z V 2 a gs φ (14)
Also, using [START_REF] Barsky | Geometric continuity of parametric curves[END_REF] and the Pythagorean identity we obtain

s φ = 1 - 1 n 2 z (15)
Then ( 14) is expressed completely in terms of the load factor as

R = V 2 a g n 2 z -1 (16)
Arriving to an expression relating the load factor and speed of an aircraft with the horizontal path of radius R that corresponds to these variables. On the other hand, if a pitch up motion is assumed, the considerations of φ = p = 0, u = V a and θ = 0 • for a maximum load factor, yields to

n z = 1 + qV a g (17)
which can be expressed as

q = (n z -1) g V a = V a R ( 18
)
where R is the radius of vertical curvature. In this manner, as the curvature is the inverse of the radius, a load factor for independent lateral or vertical maneuvers performed by the aircraft can be obtained. In this way, following the load factor limits, and depending on the speed of the aircraft, a maximum radius path can be computed, and in case that the radius of the path needs to be increased to achieve more accuracy in the deviation of the desired path from the control points, speed of the aircraft can be decreased to attend the issue.

VI. CONCLUSIONS In this paper, a formal approach for the design of 4D smooth trajectories for commercial aircraft has been developed. The proposed approach is based on the use of a particular family of curves, the Bezier curves, which are built from control points. These control points are shown to be related with path constraints issues as well as flyability constraints. The adoption of these trajectories should ease air traffic management in congested areas, as well as improve the performance of on-board guidance systems.

Fig. 1 .

 1 Fig. 1. Bezier curve completed with starting and ending straight lines.

Fig. 3 .

 3 Fig. 3. Curvature of the G 1 and G 2 continuity paths.

Fig. 4 .

 4 Fig. 4. 100m deviation optimized path for a triplet of control points using an auxiliary control point (Q3).

Fig. 5 .

 5 Fig. 5. 100m deviation optimized path.

TABLE 1 CONTROL

 1

			POINTS.
		X(m)	Y(m)	Z(m)
	1	0	0	10,000
	2 120,843	16,983	9,300
	3 210,332 -14,779	9,000
	4 272,744	-759	8,200
	5 388,920 -11,130	9,500
	6 478,501	12,964	9,800

Fig. 2. Generated path showing arc lengths and joints of Bezier curves.

TABLE 2 TIMES

 2 AND ARC LENGTHS.

	Arc length(m)		Time(s)
			t0	0
	l1	61,016	t1	305.1
	l2	107,536	t2	842.8
	l3	78,523	t3 1,235.4
	l4	89,990	t4 1,685.3
	l5	104,206	t5 2,206.4
	l6	46,383	t6 2,438.3

TABLE 3 DISTANCE

 3 OF CURVE FROM CONTROL POINTS.

	From	P2	P3	P4	P5
	Distance(m) 92.4823 99.7277 96.0905 79.5305

https://www.skybrary.aero/index.php/Accident_and_Serious_Incident_Reports:_LOS


Jean-Paul Cerri
email: jean-paul.cerri@math.u-bordeaux.fr

Pierre Lezowski
email: pierre.lezowski@math.univ-bpclermont.fr

COMPUTATION OF EUCLIDEAN MINIMA IN TOTALLY DEFINITE QUATERNION FIELDS

Keywords: 1991 Mathematics Subject Classification. Mathematics Subject Classification 2010: 11Y40, 11R04, 13F07 quaternion algebras, norm-Euclidean orders, norm-Euclidean minimum, algorithmic number theory

We describe an algorithm that allows to compute the Euclidean minimum (for the norm form) of any order of a totally definite quaternion field over a number field K of degree strictly greater than 1. Our approach is a generalization of previous work dealing with number fields. The algorithm was practically implemented when K has degree 2.

Introduction

If K is a number field, we denote by n its degree, by Z K its ring of integers, by Z × K its unit group, and by N K/Q : K → Q the norm form. Throughout this paper F will be a totally definite quaternion field over a number field K. Let us recall the relevant definitions, the reader may refer to [START_REF] Deuring | [END_REF][START_REF] Reiner | Maximal orders[END_REF][START_REF] Vignéras | Arithmétique des algèbres de quaternions[END_REF] for a complete theory. Let F be a quaternion algebra over a number field K, i.e. a 4-dimensional algebra over K with basis (1, i, j, k) such that i 2 = a, j 2 = b and k = ij = -ji, where a, b are non-zero elements of K. This algebra is denoted by a, b K . Let w = x + yi + zj + tk ∈ F , where x, y, z, t ∈ K. We denote by w the image of w under the canonical involution of F , which is defined by w = x-yi -zj -tk, and by nrd F/K (w) = ww its reduced norm. The algebra F is a division algebra if and only if the quadratic form nrd F/K (x + yi + zj + tk) = x 2 -ay 2 -bz 2 + abt 2 represents zero on K only trivially. In this case, we say that F is a quaternion field.

In addition, we will suppose that F is totally definite, or equivalently that every infinite place of K is ramified in F . This implies in particular that K is totally real, that a and b are totally negative and that for every x ∈ F \ {0}, nrd F/K (x) is totally positive. Let us denote by N : F -→ Q ≥0 the reduced norm map defined by N = N K/Q • nrd F/K . Definition 1.1. Let Λ be an order of F . For any ξ ∈ F , we set

m Λ (ξ) = inf λ∈Λ N (ξ -λ)
and we call it the local Euclidean minimum of Λ at ξ. We define the Euclidean minimum of Λ by M (Λ) = sup ξ∈F m Λ (ξ).

Set F R = F ⊗ Q R and let N R : F R → R ≥0 be the extended norm corresponding to N . Definition 1.2. Let Λ be an order of F . For any x ∈ F R , we set

mΛ (x) = inf λ∈Λ N R (x -λ)
and we call it the local inhomogeneous minimum of Λ at x. We define the inhomogeneous minimum of Λ by

M (Λ) = sup x∈F R mΛ (x).
Let us notice that these suprema are well-defined positive real numbers, that M (Λ) ≤ M (Λ) and that for every ξ ∈ F there exists a λ ∈ Λ such that m Λ (ξ) = N (ξ -λ) (see [START_REF] Chaubert | Minimum euclidien des ordres maximaux dans les algèbres centrales à division[END_REF][START_REF] Bayer | Euclidean minima and central division algebras[END_REF]). Let us call (P ) the following property: there exists some ξ ∈ F such that m Λ (ξ) = M (Λ). In particular, if (P ) holds, then M (Λ) = M (Λ) ∈ Q. Recall that we have the following results:

• The order Λ is right-norm-Euclidean if and only if it is left-norm-Euclidean (so that we can speak of norm-Euclidean orders) and that in this case Λ is necessarily a maximal order of F ;

• If M (Λ) < 1 then Λ is norm-Euclidean; • If M (Λ) > 1 then Λ is neither right neither left-norm-Euclidean;
• If M (Λ) = 1 and if property (P ) is satisfied then Λ is neither right neither left-norm-Euclidean; • If Λ is norm-Euclidean, then every other maximal order Λ ′ of F is also norm-Euclidean and M (Λ) = M (Λ ′ ) • If K has unit rank strictly greater than 1, i.e. n > 2 if K is totally real, then (P ) holds.

Note that these results still hold when F is not totally definite. For more details, the reader can refer to [START_REF] Chaubert | Minimum euclidien des ordres maximaux dans les algèbres centrales à division[END_REF][START_REF] Bayer | Euclidean minima and central division algebras[END_REF][START_REF] Cerri | Euclidean totally definite quaternion fields over the rational field and over quadratic number fields[END_REF].

A first natural question arises: what can be said when the degree of K is 1 or 2, i.e. when K = Q or is a real quadratic field? Do we also have M (Λ) = M (Λ) or does property (P ) hold? We will answer to this question in Section 2, showing that the equality is always satisfied and that property (P ) holds when K = Q. Showing that it also holds when K is quadratic seems out of reach, as in the number field case, which is a famous conjecture due to Barnes and Swinnerton-Dyer. As in the number field case (see [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF][START_REF] Lezowski | Computation of the Euclidean Minimum of Algebraic Number Fields[END_REF]), it is also natural to ask whether it is possible to use an algorithm allowing to compute M (Λ) = M (Λ) and to check that property (P ) is satisfied, even in the conjectural case. We will see that such an algorithm is already well known when K = Q so that we have just to study the case where K has degree at least 2. Note also that, as in the number field case, such an algorithm might allow to determine the upper part of the so called inhomogeneous spectrum associated to Λ. The paper is organized as follows. In Section 2, we recall what we know when n = 1 and we establish some preliminary results allowing to prove M (Λ) = M (Λ) when n = 2, and that we will use later in the algorithm. In Section 3, we explain the ideas that will be used in the algorithm when n ≥ 2. In Section 4, we describe the algorithm itself (Algorithm 2); in Section 5, we give technical details of the algorithm, and finally, in Section 6 we give the results obtained when K is quadratic, the only case that we can treat with the computers that we used.

Preliminary results

2.1.

Φ(ξ) = σ 1 (α), . . . , σ n (α), σ 1 (β), . . . , σ n (β), σ 1 (γ), . . . , σ n (γ), σ 1 (δ), . . . , σ n (δ) ,
and in further computations we see this vector as a column vector. Then Φ(Λ) is a lattice of R 4n that will be denoted by R.

Identifying F R with R 4n , if x ∈ R 4n we put N R (x) = n l=1 nrd l (x)
where

nrd l (x) = x 2 l -σ l (a)x 2 l+n -σ l (b)x 2 l+2n + σ l (ab)x 2 l+3n . With this notation, if ξ ∈ F and x ∈ R 4n , we have • nrd l (Φ(ξ)) = σ l (nrd F/K (ξ)) for every l; • N (ξ) = N R (Φ(ξ)); • mΛ (x) = inf λ∈Λ N R (x -Φ(λ)); • m Λ (ξ) = mΛ (Φ(ξ)). The multiplicative group K \ {0} acts on R 4n in the following way. If α ∈ K \ {0} and x ∈ R 4n we put α ⋆ x = (σ 1 (α)x 1 , . . . , σ n (α)x n , . . . , σ 1 (α)x 3n+1 , . . . , σ n (α)x 4n ). Hence, if ξ ∈ F we have α ⋆ Φ(ξ) = Φ(αξ). Moreover, if ν ∈ Z × K , then ν ⋆ R = R, and if x ∈ R 4n , then N R (ν ⋆ x) = N R (x). This implies that that for every x ∈ R 4n and every ν ∈ Z × K , we have mΛ (ν ⋆ x) = mΛ (x). Remark 2.1. Since mΛ is R-periodic, we have mΛ (ν ⋆ x+ Φ(λ)) = mΛ (x) for every (x, ν, λ) ∈ R 4n × Z × K × Λ.
Recall that mΛ is not only R-periodic, but also upper semi-continuous. These two properties imply that there exists an x ∈ R 4n such that M (Λ) = mΛ (x). Such an x will be called critical.

2.2.

The case n = 1. Here K = Q, a, b ∈ Q <0 and we can write

Λ = ⊕ 4
l=1 (a l,1 + a l,2 i + a l,3 j + a l,4 k)Z, where a l,i ∈ Q for every 1 ≤ l, i ≤ 4 and where the matrix M = (a l,i ) is invertible. We also have

F = ⊕ 4 l=1 (a l,1 + a l,2 i + a l,3 j + a l,4 k)Q = Q ⊕ Qi ⊕ Qj ⊕ Qk. With our notation, R is the lattice M t Z 4 ⊂ Q 4 and Φ(F ) = M t Q 4 . Let x ∈ R 4 and let us put x = M t y. Then mΛ (x) = inf X∈Z 4 q(M t (y -X)) where q(z) = z 2 1 -az 2 2 -bz 2 3 + abz 2 4
which is equivalent to

mΛ (x) = inf X∈Z 4 SM t (y -X) 2
where • is the usual Euclidean norm and

S =     1 0 0 0 0 √ -a 0 0 0 0 √ -b 0 0 0 0 √ ab     . As a consequence M (Λ)
is the square of the covering radius of the lattice L = SM t Z 4 ⊂ R 4 and if x ∈ R 4 is a critical point, i.e. satisfies mΛ (x) = M (Λ), then OP , where O = (0, 0, 0, 0) and P = Sx = SM t y, is the covering radius of L. Let A and B be two points of L and let us put A = SM t α, B = SM t β where α, β ∈ Z 4 . Then the mediator hyperplane of A and B is the set of z = SM t r such that SM t (r -α+β 2 ), SM t (α -β) = 0, i.e. r -α+β 2 , M S 2 M t (α -β) = 0, where • , • is the usual Euclidean scalar product. Note that α+β 2 , M S 2 M t (α -β) ∈ Q 4 . Now, the point P = SM t y belongs to 4 independent such mediator hyperplanes and the y i satisfy 4 independent linear equations with rational coefficients, so that y ∈ Q 4 . From this we deduce that x = M t y ∈ Φ(F ). This implies that property (P ) is satisfied and that M (Λ) = M (Λ) ∈ Q. Let us conclude by noting that there exist well known algorithms to compute the covering radius of L (see [START_REF] Sikirić | Complexity and algorithms for computing Voronoi cells of lattices[END_REF][START_REF] Schürmann | Computational Geometry of Positive Definite Quadratic Forms[END_REF]) and then M (Λ) = M (Λ), so that the case n = 1 is completely settled.

2.3.

Computation of m Λ (ξ). From now on, we suppose that n ≥ 2. We are interested in computing m Λ (ξ) where ξ ∈ F . We will rely on it for the actual computation of M (Λ) later on. Let (ε i ) 1≤i≤n-1 be a fundamental system of units of K. We denote by L the logarithmic embedding of K \ {0} in R n defined by

L(α) = (ln |σ 1 (α)|, . . . , ln |σ n (α)|).
We know that L(Z × K ) is a lattice of the hyperplane H of R n defined by the equation 1≤i≤n x i = 0, which admits (L(ε i )) 1≤i≤n-1 as a Z-basis and that the kernel of L is {±1}. For 1 ≤ i ≤ n we set

Γ i = n-1 j=1 max |σ i (ε j )|, 1 |σ i (ε j )| .
Lemma 2.2. Units of K have the following properties.

(i) Let l ∈ {1, . . . , n}. There exists a ν ∈ Z × K such that |σ i (ν)| < 1 for every i = l. (ii) If c 1 , . . . , c n-1 are given positive real numbers, there exists a ν ∈ Z × K such that c i ≤ σ i (ν) 2 ≤ c i Γ 2
i for all i ∈ {1, . . . , n -1}.

Proof. (i) L(Z × K ) being a lattice of H, it is easy to see that H ∩ {x ∈ R n ; x l > 0 and x i < 0 for i = l}, which is a nonempty open cone of H, contains an element L(ν) of L(Z × K ). (ii) This a consequence of [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF]Lemma 3.1]. See also [START_REF] Cerri | Spectres euclidiens et inhomogènes des corps de nombres[END_REF]Lemma 2.1] for a detailed proof.

Proposition 2.3. Let k > 0. Suppose that x ∈ R 4n and X ∈ R satisfy N R (x-X) < k. Then there exist an ε ∈ Z × K and some Y ∈ R such that y = ε ⋆ x -Y satisfies N R (y) < k and 0 ≤ nrd i (y) ≤ Γ(k) for every 1 ≤ i ≤ n,
where

Γ(k) = k n-1 i=1 Γ 2 i 1 n . Proof. First, if N R (x -X) = 0 we have nrd l (x -X) = 0 for some index l. By Lemma 2.2 (i), there exists a ν ∈ Z × K such that |σ i (ν)| < 1 for every i = l. Since nrd i (ν p ⋆ x -ν p ⋆ X) = σ i (ν) 2p nrd i (x -X), for p sufficiently large we have nrd i (ν p ⋆ x -ν p ⋆ X) ≤ Γ(k) for i = l and nrd l (ν p ⋆ x -ν p ⋆ X) = 0. Taking ε = ν p and Y = ν p ⋆ X, we obtain the result announced. Now, if N R (x -X) = 0, we apply Lemma 2.2 (ii) with c i = Γ(k) Γ 2 i nrd i (x -X) for 1 ≤ i ≤ n -1. There exists a ν ∈ Z × K such that Γ(k) Γ 2 i nrd i (x -X) ≤ σ i (ν) 2 ≤ Γ(k) nrd i (x -X) for 1 ≤ i ≤ n -1. From N R (ν ⋆ x -ν ⋆ X) = N R (x -X) = σ n (ν) 2 nrd n (x -X) n-1 i=1 σ i (ν) 2 nrd i (x -X) we deduce σ n (ν) 2 nrd n (x -X) ≤ k n-1 i=1 Γ 2 i Γ(k) n-1 = Γ(k).
Taking ε = ν and Y = ν ⋆ X we have again the result announced. Now, consider the orbits of elements of R 4n under the action of

Z × K . Let x, y ∈ R 4n and ν ∈ Z × K . As x -y ∈ R ⇒ ν ⋆ x -ν ⋆ y ∈ R, the group Z × K also acts on R 4n /R by (ν, x) → ν ⋆ x, where y is the class of y ∈ R 4n modulo R. Let F be a fundamental domain of R. Identifying R 4n /R with F , if x ν is the unique element of F congruent to ν ⋆ x modulo R, we set Orb(x) = {x ν ; ν ∈ Z × K }. Remark that for all x ∈ R 4n , if z ∈ Orb(x) we have mΛ (z) = mΛ (x).
Using this new notation, we have the following essential result.

Theorem 2.1. There exists a finite set S ⊂ R such that for every x ∈ R 4n , we have

mΛ (x) = inf z∈Orb(x) min Z∈S N R (z -Z) .
Proof. Let k ′ be any positive real number satisfying

k ′ > M (Λ) and let ǫ > 0 such that M (Λ) + ǫ < k ′ . Let S ′ = {t ∈ R 4n ; nrd i (t) ≤ Γ(k ′ ) for every 1 ≤ i ≤ n}.
Since for every i,

nrd i (t) = t 2 i -σ i (a)t 2 i+n -σ i (b)t 2 i+2n + σ i (ab)t 2 i+3n , where σ i (a), σ i (b) < 0, we see that S ′ is a bounded subset of R 4n . Now, let us put S = (F + S ′ ) ∩ R. As S ′ + F is also bounded, the set S is a finite subset of R. Let x ∈ R 4n . There exists some X ∈ R such that N R (x -X) < mΛ (x) + ǫ ≤ M (Λ) + ǫ < k ′ . Applying Proposition 2.3 with k = mΛ (x) + ǫ, we see that there exists a ν ∈ Z × K and some Y ∈ R such that, if y = ν ⋆ x -Y , then N R (y) < mΛ (x) + ǫ < k ′ and nrd i (y) ≤ Γ( mΛ (x) + ǫ) ≤ Γ(k ′ ) (because Γ is an increasing function). Let z ∈ F such that z -y ∈ R. Then, z ∈ Orb(x) and y = z -Z, where Z ∈ S and N R (z -Z) < mΛ (x) + ǫ. Since this is true for any ǫ > 0 we have mΛ (x) ≥ inf z∈Orb(x) min Z∈S N R (z -Z) .
On the other hand, for every z ∈ Orb(x) and any Z ∈ R, we have mΛ (x) = mΛ (z) ≤ N R (z -Z), which leads to the equality.

Remark 2.4. This result shows that if x ∈ Φ(F ), we can compute mΛ (x) in a finite number of steps. In fact, Orb(x) is finite if and only x ∈ Φ(F ). Indeed, if x ∈ Φ(F ) we can write x = X/d where X ∈ R and d ∈ Z > 0, and Orb(x) can be identified with a subset of R/dR which is finite. Conversely, if Orb(x) is finite and if ν is a nontorsion unit of K, considering the sequence (x ν p ) p≥0 , we see that there exists k > l ≥ 0 such that x ν k = x ν l . This implies that there exists an

X = Φ(λ) ∈ R, with λ ∈ Λ, such that (ν k -ν l ) ⋆ x = X which implies x = Φ(λ/(ν k -ν l )) ∈ Φ(F ).
As in the number field case, this allows to establish an algorithm to compute mΛ (x). We could use the formula for S in the proof of Theorem 2.1 but we do not know M (Λ), so we cannot compute S. However, we can proceed as follows. Let us take some k ′ > 0 and set S k ′ = S as defined in the previous proof. Let us put (1)

M k ′ = min z∈Orb(x) min Z∈S k ′ N R (z -Z) .
Then

M k ′ ≤ k ′ ⇒ mΛ (x) = M k ′ .
Of course, mΛ (x) ≤ M k ′ . Suppose that mΛ (x) < M k ′ so that there exists some

X ∈ R such that N R (x -X) < M k ′ . Then by Proposition 2.3 there exists a ν ∈ Z × K and some Y ∈ R such that if y = ν ⋆ x -Y , then N R (y) < M k ′ and nrd i (y) ≤ Γ(M k ′ ) ≤ Γ(k ′ ) for every i. This contradicts the definition of M k ′ .
Finally, we have the following algorithm to compute mΛ (x).

Algorithm 1: Computation of mΛ (x) when x ∈ Φ(F ) Input: x ∈ Φ(F ). Output: mΛ (x). 1 Computation of Orb(x). 2 Computation of k ′′ = M k ′ for a k ′ > 0. 3 if k ′′ ≤ k ′ then 4 mΛ (x) = k ′′ 5 else 6 Computation of k = M k ′′ . 7 mΛ (x) = k.
The correctness of the algorithm follows from the fact that k

′ -→ M k ′ is a decreasing function: if k ′′ > k ′ then M k ′′ ≤ M k ′ = k ′′ and mΛ (x) = M k ′′ . 2.4. The case n = 2.
Here we use Theorem 2.1 to show that when n = 2, as in every other case, we have

M (Λ) = M (Λ). Theorem 2.2. If n = 2, then M (Λ) = M (Λ).
Proof. We follow the idea of the proof given by Barnes and Swinnerton-Dyer in the real quadratic number field case [START_REF] Barnes | The inhomogeneous minima of binary quadratic forms (II)[END_REF].

Let K = Q( √ d), where d > 1 is a squarefree integer. Denote by σ the nontrivial embedding of K into R, i.e. the Q-automorphism of K such that σ( √ d) = - √ d.
Let ν be a fundamental unit of K. We can suppose ν > 1 so that |σ(ν)| < 1. Let ǫ > 0. Keep the notation of the proof of Theorem 2.1. The map N R being continuous, let ǫ ′ > 0 such that for every x, y

∈ F + S + {z ∈ R 4n ; |z i | < 1 for all i}, which is a bounded set, if |x i -y i | < ǫ ′ for every i, then |N R (x) -N R (y)| < ǫ. Choose ǫ ′ < 1. Let x ∈ R 8 such that mΛ (x) = M (Λ). If x ∈ Φ(F )
we are done, so that we can suppose x ∈ Φ(F ). The set {x ν p ; p ≥ 0} ⊂ F is infinite and bounded , and there exist

n 2 > n 1 ≥ 0 such that for every i, |(x ν n 2 ) i -(x ν n 1 ) i | < ǫ ′ /2. Without loss of generality we can suppose ν n2-n1 > 2. We have ν n2-n1 ⋆ x ν n 1 = x ν n 2 + X where X ∈ R. Let us put X = Φ(λ) with λ ∈ Λ. Set ξ = λ ν n2-n1 -1 ∈ F and α = Φ(ξ) ∈ Φ(F ). Since ν n2-n1 ⋆ α = α mod R, we have Orb(α) = n2-n1-1 k=0 {α ν k , α -ν k } ⊂ F .
It is easy to see that

α i =      (x ν n 1 ) i + (x ν n 1 ) i -(x ν n 2 ) i ν n2-n1 -1 for i odd, (x ν n 1 ) i + (x ν n 1 ) i -(x ν n 2 ) i σ(ν) n2-n1 -1 for i even. Since |(x ν n 1 ) i -(x ν n 2 ) i | < ǫ ′ /2 and ν n2-n1 > 2, this implies that |α i -(x ν n 1 ) i | < ǫ ′ ν n2-n1 for i odd, and |α i -(x ν n 1 ) i | < ǫ ′ for i even. Consequently for every 0 ≤ k < n 2 -n 1 , we have (2) |(ν k ⋆ α) i -(ν k ⋆ x ν n 1 ) i | < ǫ ′ for all i. Now, for every 0 ≤ k < n 2 -n 1 , ν k ⋆ α = α ν k + X k where X k ∈ R. Let X ∈ S.
Then for every i, we have

|(α ν k -X) i -(ν k ⋆ x ν n 1 -X k -X) i | < ǫ ′ . But α ν k -X ∈ F + S and since ǫ ′ < 1, we have ν k ⋆ x ν n 1 -X k -X ∈ F + S + {z ∈ R 4n ; |z i | < 1 for all i}. By choice of ǫ ′ this implies that for every 0 ≤ k < n 2 -n 1 , (3) |N R (α ν k -X) -N R (ν k ⋆ x ν n 1 -X k -X)| < ǫ.
Now inequality [START_REF] Bayer | Euclidean minima and central division algebras[END_REF] shows that for every 0

≤ k < n 2 -n 1 , |(-ν k ⋆ α) i -(-ν k ⋆ x ν n 1 ) i | < ǫ ′ for all i. Putting -ν k ⋆ α = α -ν k + Y k where Y k ∈ R for 0 ≤ k < n 2 -n 1 ,
we obtain in the same way that for every x ∈ S and every 0

≤ k < n 2 -n 1 , (4) |N R (α -ν k -X) -N R (-ν k ⋆ x ν n 1 -Y k -X)| < ǫ.
Thus, by ( 3) and ( 4), for every y ∈ Orb(α) and every X ∈ S,

N R (y -X) > N R (z -Z) -ǫ
for some z ∈ Orb(x ν n 1 ) and Z ∈ R. We deduce from Theorem 2.1 that

mΛ (α) > mΛ (x ν n 1 ) -ǫ = mΛ (x) -ǫ = M (Λ) -ǫ.
Since for every ǫ > 0, there exists some α ∈ Φ(F ) such that mΛ (α) > M (Λ) -ǫ, we have necessarily M (Λ) = M (Λ).

Remark 2.5. The second part of our proof does not follow [START_REF] Barnes | The inhomogeneous minima of binary quadratic forms (II)[END_REF] because the third "equation" of [1, p. 313] is incorrect.

Corollary 2.6. In all cases, we have M (Λ) = M (Λ).

Proof. By [START_REF] Bayer | Euclidean minima and central division algebras[END_REF], the equality holds when n > 2, and by Subsection 2.2, it also holds when n = 1. Of course, in both cases we have a better result since property (P ) holds.

Theoretic argumentation

3.1. Overview of the strategy. Now, it is time to set out the ideas which are behind the algorithm used to compute M (Λ) = M (Λ) when n ≥ 2. The strategy is the same as in the number field case.

To simplify things we assume that we have an idea of the exact value of M (Λ). We shall see later how one can find a good candidate for the Euclidean minimum. From now on, we denote by k our guess of M (Λ).

In fact, instead of proving the equality M (Λ) = k, we shall establish the stronger and more precise result:

M (Λ) ≤ k and there exists a ξ ∈ F such that m Λ (ξ) = k.
It will clearly follow that M (Λ) = M (Λ) = k. Moreover, we shall try to find all the critical points which belong to F . Since mΛ is R-periodic, it is sufficient to work on a fundamental domain F of R, i.e. to prove that for all x ∈ F , mΛ (x) ≤ k, and to find all the ξ ∈ F such that Φ(ξ) ∈ F and m Λ (ξ) = k (every solution to m Λ (ξ) = k will be of this form modulo Λ).

Let k ′ be a positive number smaller than k. In practice one takes k ′ = k -ǫ where ǫ is a small positive number. Let us consider a finite family of elements of R, say X , and the regions of R 4n centered in the X of X and defined by the inequations N R (x -X) ≤ k ′ . Every element x of the subset H of F covered by these regions satisfies mΛ (x) ≤ k ′ < k, and since k ′ is supposed smaller than M (Λ), "holes" appear in the covering of F by these regions. These holes contain the potentially critical points of F .

The main idea is then to analyze the action of the unit group Z × K on uncovered subsets of F , in the following way. Let T be a hole of F , and ε a non-torsion unit of K (ε = ±1). We look at the possible intersections of ε ⋆ T with holes of F modulo R. If ε ⋆ T does not intersect any hole of F modulo R, we know by Remark 2.1 that for every x of T , we have mΛ (x) ≤ k ′ , so that T can be eliminated as a subset of F potentially containing a critical point. The interesting case is when the intersection is nonempty.

Remark 3.1. Here we have expressed things in terms of holes. In what follows, we consider "easy" regions larger than holes. For instance, in the algorithm, holes are replaced by regions composed of small parallelotopes. All what we need is to have a partition of F in a covered region and in regions potentially containing critical points. Then we check that these regions have an exploitable behaviour under the action of Z × K . 3.2. Theoretic arguments. As in the previous subsection, we consider k ′ > 0 and and a subset H of R 4n which satisfies mΛ (x) ≤ k ′ for all x ∈ H. We consider also a unit ε = ±1.

The cyclic case.

Let us first study the cyclic situation.

Theorem 3.1. Let T 0 ,. . . , T j-1 be bounded subsets of R 4n (j ≥ 1). Assume that for all l there is an Υ l ∈ Λ such that

(5) (ε ⋆ T l -Φ(Υ l )) \ H ⊂ T l+1 ,
where the indices in T r are to be read modulo j. Assume also that there is an x in T 0 which satisfies mΛ (x) > k ′ and define Ω ∈ Λ by

Ω = ε j-1 Υ 0 + ε j-2 Υ 1 + . . . . + εΥ j-2 + Υ j-1 .
Consider the sequence defined by y 0 = x and y p+1 = ε j ⋆ y p -Φ(Ω) for all p ≥ 0.

Then, if we put

ξ = Ω ε j -1 ∈ F and t = Φ(ξ),
we have i) For all i ∈ {1, . . . , n} such that |σ i (ε)| > 1 and for all p ≥ 0 and all k ∈ {0, 1, 2, 3}, (y p

) i+kn = t i+kn . ii) The sequence (y p ) p≥0 converges to t. iii) k ′ < mΛ (x) ≤ mΛ (t). iv) If x ∈ Φ(F ) then x = t. Proof. Put H ′ = {x ∈ R 4n ; mΛ (x) ≤ k ′ } so that H ⊂ H ′ . First of all, let us prove that (6) ε j ⋆ T 0 -Φ(Ω) \ H ′ ⊂ T 0 .
Put z = ε j ⋆ z 0 -Φ(Ω) where z 0 ∈ T 0 and suppose z ∈ H ′ . Let us define z 1 , z 2 , . . . , z j by the induction formula z p+1 = ε ⋆ z p -Φ(Υ p ) for 0 ≤ p < j.

It is easy to see that we have z j = z. By Remark 2.1 we can write

mΛ (z) = mΛ (z j ) = mΛ (z j-1 ) = . . . = mΛ (z 0 ) > k ′ .
Thus for all p ∈ {0, . . . , j} we have z p ∈ H ′ , which implies z p ∈ H, and by successive applications of (5) we get z p ∈ T p for all p ∈ {0, . . . , j -1}, and finally z = z j ∈ T 0 , so that we have [START_REF] Cerri | Euclidean totally definite quaternion fields over the rational field and over quadratic number fields[END_REF]. Now, consider the sequence (y p ) p≥0 . By Remark 2.1 we see by induction that for all p ≥ 0

(7) mΛ (y p ) = mΛ (x) > k ′
so that for all p ≥ 0, y p ∈ H ′ . Then, as y 0 = x ∈ T 0 , using [START_REF] Cerri | Euclidean totally definite quaternion fields over the rational field and over quadratic number fields[END_REF] we easily establish by induction that y p ∈ T 0 for all p ≥ 0. Thus, as T 0 was assumed to be bounded, the sequence (y p -t) p≥0 is bounded.

But, by the definition of t and the induction formula which defines (y p ) p≥0 , we have y p -t = ε jp ⋆ (x -t) for all p ≥ 0, so that

(8) |(y p ) i -t i | = |σ i mod n (ε)| jp |x i -t i |
for all i ∈ {1, 2, . . . , 4n} and for all p ≥ 0.

Let i ∈ {1, . . . , n}. If |σ i (ε)| > 1, since the sequence (|(y p ) i+kn -t i+kn |) p≥0 is bounded for every k ∈ {0, 1, 2,
3}, we must have x i+kn -t i+kn = 0, and then by ( 8) we obtain [START_REF] Deuring | [END_REF] (y p ) i+kn = t i+kn for all k ∈ {0, 1, 2, 3} and for all p ≥ 0.

This is i).

Moreover if |σ i (ε)| < 1, then [START_REF] Chaubert | Minimum euclidien des ordres maximaux dans les algèbres centrales à division[END_REF] shows that (10) lim p→+∞ (y p ) i+kn = t i+kn for all k ∈ {0, 1, 2, 3}.

Since |σ i (ε)| = 1 (otherwise ε = ±1, which is excluded by hypothesis), ( 9) and ( 10) yield lim p→+∞ y p = t. This is ii). Finally, since mΛ is upper semi-continuous, by [START_REF] Cerri | Tables of Euclidean minima of orders in totally definite quaternion fields over quadratic fields[END_REF], we obtain:

k ′ < mΛ (x) = lim sup p→+∞ mΛ (y p ) ≤ mΛ (t),
which gives iii). Now, assume that x ∈ Φ(F ). Since we cannot have |σ i (ε)| ≤ 1 for every i, there exists an i ∈ {1, . . . , n} such that |σ i (ε)| > 1, and by i) (with p = 0) we have [START_REF] Lezowski | Computation of the Euclidean Minimum of Algebraic Number Fields[END_REF] x i+kn = t i+kn for every k ∈ {0, 1, 2, 3}.

But x and t are both in Φ(F ), and if we write x = α 0 + α 1 i + α 2 j + α 3 k, t = β 0 + β 1 i + β 2 j + β 3 k with α j , β j ∈ K for every j, by [START_REF] Lezowski | Computation of the Euclidean Minimum of Algebraic Number Fields[END_REF] we have σ i (α j ) = σ i (β j ) for every j. By injectivity of σ i , this leads to α j = β j for every j, so that x = t.

Remark 3.2. Obviously the same property holds for T 1 , T 2 , . . . , T j-1 , the single difference is the formula for Ω, in which indices must be trivially permuted. More precisely, for r ∈ {0, . . . , j -1}, if we put t r = Φ(ξ r ) with

ξ r = Ω r ε j -1 ,
and

Ω r = ε j-1 Υ r + ε j-2 Υ r+1 + . . . + Υ r+j-1 ,
where the indices are still to be read modulo j, we have the same property as in Theorem 3.1 for T r (with t r instead of t). Moreover t 0 = t and we have the cyclic law:

t r+1 = Φ(ε) • t r -Φ(Υ r
) for all r ∈ {0, . . . , j -1}. In particular, all the t r are in Orb(t).

3.2.2. Generalization. Even if Theorem 3.1 allows one to treat some situations, it is not sufficient, in the form seen above, to cover all the cases that one meets in practice. A generalization of the previous situation is the following one.

Let T i (0 ≤ i ≤ s -1) be distinct bounded sets of R 4n , and T = {T 0 , . . . , T s-1 }. Assume that for all T i in T there exists an X i ∈ R and s i integers n i,1 ,. . . ,n i,si (s i > 0) such that [START_REF]PARI/GP version 2.9[END_REF] (ε

⋆ T i -X i ) \ H ⊂ 1≤k≤si T n i,k .
To simplify notation we shall consider the T i as the vertices of a directed graph (from now on digraph) G and represent ( 12) by s i directed edges (from now on arcs) whose tail is T i and whose respective heads the

T n i,k (1 ≤ k ≤ s i ).
Of course, such an arc can be a loop. We shall write

T i → T n i,k (X i ) or T i → T n i,k , if it is not necessary to precise X i .
Example 3.3. Theorem 3.1 corresponds to the digraph

G 1 : T 0 → T 1 (Φ(Υ 0 )) → . . . → T j-1 (Φ(Υ j-2 )) → T 0 (Φ(Υ j-1 )).
To describe paths of G we shall use the notation

T ′ 1 → T ′ 2 → . . . → T ′ k .
The digraph G has the following properties: if T and T ′ are vertices of G, there is at most one arc whose tail is T and whose head is T ′ , and every vertex of G has a positive outvalency. Obviously the last property implies that G contains circular paths (or circuits). Consequently, the set C of simple circuits of G (paths of the form

T ′ 0 → . . . → T ′ k → T ′ 0 ,
where k ≥ 0 and all the T ′ i are distinct) is nonempty (take a circuit of minimal length) and is finite (their length cannot exceed s). Each element c of C of length j is of the form of the circular path met in Theorem 3.1 (and seen above in G 1 ),

T ′ 0 → T ′ 1 (X ′ 0 ) . . . → T ′ j-1 (X ′ j-2 ) → T ′ 0 (X ′ j-1 ) with X ′ i = Φ(Υ i ).
It defines, in a unique way, j points of Φ(F ), t 0 , . . . , t j-1 by the formulae of Remark 3.2.

In this context, we say that t 0 , . . . , t j-1 are associated to c (implicitly t i corresponds to T ′ i ). The t i are in the same orbit modulo R and satisfy mΛ (t 0 ) = . . . = mΛ (t j-1 ). Let us denote this rational number by m Λ (c) and put

m Λ (G) = max c∈C m Λ (c).
Moreover, let us denote by E the set of all points of Φ(F ) associated to the elements of C. The set E is finite and we also have

m Λ (G) = max t∈E mΛ (t).

Finally let us put

E ′ = {t ∈ E such that mΛ (t) = m Λ (G)}.
An infinite path of G is an infinite sequence of arcs of G, (A i ) i≥0 such that the head of A i is the tail of A i+1 . If A i is defined by T ′ i → T ′ i+1 , we shall denote the path by (T ′ i ) i≥0 . Such a path is not simple, but can have a periodicity property. An infinite path (T ′ i ) i≥0 is said to be ultimately periodic if there exist integers r ≥ 0 and p ≥ 1 such that [START_REF] Reiner | Maximal orders[END_REF] for all i ≥ r, T ′ i+p = T ′ i .

Let (T ′ i ) i≥0 be an ultimately periodic infinite path. Let P be the set of p ≥ 1 such that there exists an r which satisfies [START_REF] Reiner | Maximal orders[END_REF]. Then P is nonempty and we can define ρ = min P ≥ 1.

Then there exists an r ρ such that ∀i ≥ r ρ , T ′ i+ρ = T ′ i . The integer ρ will be called the period length of (T ′ i ) i≥0 and every circuit T ′ i → . . . → T ′ i+ρ , where i ≥ r ρ , will be called a period of (T ′ i ) i≥0 Definition 3.4. We shall say that G is convenient if every infinite path of G is ultimately periodic.

Convenient digraphs have the following properties. Assume that G is convenient. Then every circuit is a power of a simple circuit. It can be shown that this condition implies that G is convenient. Another characterisation of convenient digraphs could be the following one: two distinct simple circuits have no common vertex. Assume that G is convenient and let P = (T ′ i ) i≥0 be an infinite path of G. Then, every period of P is a simple circuit (see [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF] for more details).

Example 3.5. G 1 is convenient. Figure 1 gives examples of convenient and not convenient digraphs.

T 2 T 3 T 4 T 5 T 6 T 1 X 2 X 5 X 3 X 4 X 6 X 7 X 1 (a) G1 is convenient T 2 T 3 T 4 T 5 X 1 X 2 X 3 X 4 X 5 X 6 (b) G2 is not convenient (X2 = X3) T 1 T 2 T 3 T 4 T 5 T 6 X 1 X 1 X 1 X 1 X 2 X 3 X 4 X 5 X 6 (c) G3 is not convenient T 1 T 2 T 3 T 4 T 5 T 6 T 7 X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 9 X 10 X 8 (d) G4 is convenient Figure 1. Some digraphs.
Now, we can establish the theorem which will allow us to treat (almost) all the situations. Theorem 3.2. Assume that G is convenient and that there exist a T ∈ T and an

x ∈ T such that mΛ (x) > k ′ . Then i) k ′ < mΛ (x) ≤ m Λ (G). ii) If x ∈ Φ(F ), there exists a t ∈ E such that x ≡ t mod R. iii) If x ∈ Φ(F ) is critical, there exists a t ∈ E ′ such that x ≡ t mod R.
Proof. The proof is exactly the same as the proof of Theorem 4.5 of [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF] with some changes of notation. Let us rephrase it. Put x 0 = x and T ′ 0 = T . By [START_REF]PARI/GP version 2.9[END_REF] we know that there exists X ′ 0 ∈ R and s ′ 0 elements of T , denoted by

T n ′ 0,k (1 ≤ k ≤ s ′ 0 ) such that (ε ⋆ T ′ 0 -X ′ 0 ) \ H ⊂ 1≤k≤s ′ 0 T n ′ 0,k . Set x 1 = ε ⋆ x 0 -X ′ 0 . Since mΛ (x 1 ) = mΛ (x 0 ) > k ′ , we have x 1 ∈ (ε ⋆ T ′ 0 -X ′ 0 ) \ H, and necessarily, there is an i ∈ {1, . . . , s ′ 0 } such that x 1 ∈ T n ′ 0,i . We put T ′ 1 = T ′ n ′ 0,i
, and we continue with

x 2 = ε ⋆ x 1 -X ′ 1 where X ′ 1 is the element of R associated to T ′
1 by [START_REF]PARI/GP version 2.9[END_REF]. We see that we can construct by induction a sequence (x i ) i≥0 and an infinite path (T ′ i ) i≥0 which satisfy:

x 0 = x, for all i ≥ 0, x i+1 = ε ⋆ x i -X ′ i where X ′ i ∈ R and (14)
for all i ≥ 0, x i ∈ T ′ i . Moreover, by Remark 2.1, we have mΛ (x i ) = mΛ (x) > k ′ for all i. G being convenient, the infinite path (T ′ i ) i≥0 is ultimately periodic. We denote its period length ρ and we consider one of its periods c, described by

T ′ r → . . . → T ′ r+ρ = T ′ r , which is a simple circuit. Define T ′′ s = {x r+s+iρ ; i ∈ N}, for 0 ≤ s ≤ ρ -1.
By [START_REF] Schürmann | Computational Geometry of Positive Definite Quadratic Forms[END_REF] we have ∀s ∈ {0, . . . , ρ -1}, T ′′ s ⊂ T ′ r+s . This implies that the T ′′ s are bounded. Moreover, by construction, for all s there exists

Υ s ∈ Λ (in fact Φ -1 (X ′ r+s )) such that ε ⋆ T ′′ s -Φ(Υ s ) \ H = ε ⋆ T ′′ s -Φ(Υ s ) ⊂ T ′′ s+1
where the indices are to be read modulo ρ. Putting y = x r ∈ T ′′ 0 which satisfies mΛ (y) > k ′ , we see that we are exactly under the hypotheses of Theorem 3.1 (with y instead of x, T ′′ i instead of T i and ρ instead of j). This theorem defines ρ rational points t i associated to the simple circuit c. By definition of m Λ (c), and by Theorem 3.1.iii), we obtain

k ′ < mΛ (x) = mΛ (x r ) ≤ m Λ (c),
and by definition of m Λ (G) we have i).

Assume now that x ∈ Φ(F ) so that, by induction, x r ∈ Φ(F ). By Theorem 3.1.iv) we have x r = t 0 , and thus x r -t 0 ∈ R. By the induction formula of the definition of (x i ) and the formulae of Remark 3.2, we see that

for all k ∈ {0, . . . , r}, ε ⋆ (x r-k -t -k ) ∈ R,
where the index in t -k is taken modulo ρ. Finally, x = x 0 ≡ t -r mod R which is an element of E by definition of E. This proves ii).

Assume now that x is critical so that we have mΛ (x) = M (Λ). From the definitions, we can write mΛ (x) ≥ m Λ (G) and by i) we obtain mΛ (x) = m Λ (G) so that mΛ (t -r ) = m Λ (G). Since t -r ∈ E, we find t -r ∈ E ′ . This proves iii).

4. The algorithm, theoretic aspect 4.1. General strategy. Let Λ be an order of F . Since Λ is a torsion-free module over Z which is a PID, Λ admits a Z-basis whose cardinality is 4n. Suppose that we know such a basis (e 1 , e 2 , . . . , e 4n ). For every 1 ≤ l ≤ 4n, let us write e l = β l,i + β l,2 i + β l,3 j + β l,4 k where β l,m ∈ K for every m. Now, consider the matrix M ∈ M 4n (R) defined by

M =                           σ 1 (β 1,1 ) σ 1 (β 2,1 ) . . . σ 1 (β 4n,1 ) σ 2 (β 1,1 ) σ 2 (β 2,1 ) . . . σ 2 (β 4n,1 ) . . . . . . . . . σ n (β 1,1 ) σ n (β 2,1 ) . . . σ n (β 4n,1 ) σ 1 (β 1,2 ) σ 1 (β 2,2 ) . . . σ 1 (β 4n,2 ) . . . . . . . . . σ n (β 1,2 ) σ n (β 2,2 ) . . . σ n (β 4n,2 ) σ 1 (β 1,3 ) σ 1 (β 2,3 ) . . . σ 1 (β 4n,3 ) . . . . . . . . . σ n (β 1,3 ) σ n (β 2,3 ) . . . σ n (β 4n,3 ) σ 1 (β 1,4 ) σ 1 (β 2,4 ) . . . σ 1 (β 4n,4 ) . . . . . . . . . σ n (β 1,4 ) σ n (β 2,4 ) . . . σ n (β 4n,4 )                           .
Then it is not difficult to see that we have

Φ(F ) = M •Q 4n , R = Φ(Λ) = M •Z 4n , F = M • [0, 1) 4n , ∆ = M • (Q ∩ [0, 1)) 4n and F = ∆ = M • [0, 1] 4n
. Now, as in the totally real number field case (see [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF]) we consider a cutting-covering of F = ∆ using parallelotopes whose faces are orthogonal to the canonical axes of R 4n . These parallelotopes P are of the form ( 15)

P = {(u l ) 1≤l≤4n ∈ R 4n ; |u l -C l | ≤ h l },
where C = (C l ) 1≤l≤4n is the center of the parallelotope and 0 < h l for every l.

The way we obtain a cutting-covering of F is similar to the way we proceed in the totally real number field case, so we refer to [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF] for details. The structure of the algorithm is also the same. Let us recall its general strategy and explain where differences occur. Assume that we have an idea of M (Λ) denoted k. Suppose that we have at our disposal a set X of elements of R, and let us take a small ǫ > 0. Definition 4.1. A subset of R 4n will be said to be absorbed by X ∈ X , if it is contained in the region defined by the inequality N R (x -X) ≤ k -ǫ.

The computations are organized as described in Algorithm 2. 4.2. The absorption test. Let P be a parallelotope defined as in [START_REF] Vignéras | Arithmétique des algèbres de quaternions[END_REF] and X be some element of R. Since for t ∈ R 4n we have

N R (t) = n i=1 (t 2 i -σ i (a)t 2 i+n -σ i (b)t 2 i+2n + σ i (ab)t 2 i+3n ),
where for every i, σ ( a), σ i (b) < 0, we see that P is absorbed by X if (16)

n i=1 A i (P, X) ≤ k -ǫ.
where

A i (P, X) = |C i -X i | + h i 2 -σ i (a) |C i+n -X i+n | + h i+n 2 -σ i (b) |C i+2n -X i+2n | + h i+2n 2 +σ i (ab) |C i+3n -X i+3n | + h i+3n 2 .
where for every i,

ν i ∈ Z ∩ [⌊α i ⌋, ⌊β i ⌋], α i = 4n j=1 m ′ i,j (C ′ j -δ i,j h ′ j ), β i = 4n j=1 m ′ i,j (C ′ j + δ i,j h ′ j
) and δ i,j = 1 or -1 according to whether m ′ i,j > 0 or not. Remark that X 0 satisfies (17) and denote by X 0 , X 1 , . . . , X g (g ≥ 0) all the elements of R satisfying (17). As in [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF], if for all j (0 ≤ j ≤ g) we have for all p ∈ {1, . . . , N }, there exists i p ∈ {1, . . . n} and k ∈ {0, . . . , 3} such that |(C ′′ -X j -D p ) ip+kn | > 1 + |σ ip (ε)| h ip+kn , then P can be eliminated from the list of problematic parallelotopes. 4.4. Termination without failure. It can be shown, using arguments similar to those that are used in [START_REF] Bayer | Euclidean minima and central division algebras[END_REF][START_REF] Chaubert | Minimum euclidien des ordres maximaux dans les algèbres centrales à division[END_REF][START_REF] Cerri | Spectres euclidiens et inhomogènes des corps de nombres[END_REF] (see in particular [START_REF] Cerri | Spectres euclidiens et inhomogènes des corps de nombres[END_REF]Proposition 4.25]), that the algorithm terminates as soon as n > 2 for a good choice of k and a sufficiently thin initial cutting-covering. But, in the case n = 2, which is the only case that we can treat with our computers, nothing can be said a priori.

5.

The algorithm, technical aspect 5.1. Symmetry and fundamental domain. For any x ∈ R 4n , we have mΛ (x) = mΛ (-x). Therefore, we can restrict our calculation to one half of F in Step 3 of Algorithm 2. Nevertheless, the symmetry should be kept in mind for the remainder of the algorithm (e.g. for Step 4, we should take it into account for G).

Some problematic parallelotopes can lie on the boundary of the fundamental domain, which may add unnecessary vertices to the digraph, making it non convenient. To tackle this issue, we change the fundamental domain by translating the problematic parallelotopes or the fundamental domain itself (which is equivalent anyway). For more details, see [4, Section 5.7.2]. 5.2. Initial value of k. Algorithm 2 requires an initial value of k in Step 2. For the algorithm to be successful, we need k ≤ M (Λ). If every parallelotope is eliminated by the absorption of integers or the action of units, then M (Λ) < k and we can start again with smaller k.

We observe that the points x ∈ F such that M (Λ) = m Λ (x) have a "small" orbit under the action of units. Therefore, if we denote by ε a unit of K which is not a root of unity, then ε ℓ x -x ∈ Λ for some "small" ℓ ∈ Z >0 . In other words, there exist some a ∈ Λ and some "small" ℓ ∈ Z >0 such that x = a ε ℓ -1 .

We may study such points and apply Algorithm 1 to find their minimum m Λ (x). This can provide a lower bound on M (Λ). This heuristic can turn out to be non applicable when N K/Q (ε ℓ -1) is large. Instead, we can simply consider points

x = a b ,
where b is a factor of ε ℓ -1 of small norm.

5.3.

Analysis of the digraph. In Step 6 of Algorithm 2, we try to apply Theorem 3.2, that is to say to obtain a convenient graph. This digraph may be not convenient, but we can try to simplify it to make it convenient, as in [11, Section 3.2.2]. The basic idea is to merge parallelotopes which have the same translation vector. For example, in Figure 1c, we can merge T 1 , T 2 , and T 3 to make the digraph convenient.

Execution of Algorithm 1.

Once we have a convenient digraph, it remains to apply Algorithm 1 to a finite set of points. This algorithm is described with elements of R 4n , which would suggest using floating point numbers, but it actually deals with elements of F . Consequently, the orbit is computed as a subset of the quaternion field F in Step 1.

In practice, we compute the orbit of x = Φ(ξ) ∈ Φ(F ), where ξ ∈ F as follows. Let {ε i , 1 ≤ i < n} be a fundamental system of units of K. For fixed 1 ≤ i < n and for increasing k ≥ 1, we compute the elements

ε k i • ξ reduced modulo Λ, which is denoted by ε k i • ξ red , until ε k i • ξ red = ξ
red , which happens for some finite m i ≥ 1.

Then the orbit of x is the image by Φ of the set

   ± n-1 i=1 ε ki i • ξ red red , 0 ≤ k j < m j for any 1 ≤ j < n    ,
where repetitions are discarded. As for Steps 2 and 6, to compute M k ′ or M k , we deal with quantities N R (z -Z) (see [START_REF] Barnes | The inhomogeneous minima of binary quadratic forms (II)[END_REF]) which are in fact rational numbers.Therefore, during the actual execution of the algorithm we can compute exactly these numbers and the output of Algorithm 2 is a rational number (whenever its execution is successful).

Results

6.

1. An example. Let us introduce this section with an example. We know by [START_REF] Bayer | Euclidean minima and central division algebras[END_REF]

that if K = Q( √ 2), F = -1, -1 K and Λ is a maximal order of F , then M (Λ) ≤ 1 2 .
Let us consider the maximal order

Λ = Z K ⊕ Z K √ 2(1 + i) 2 ⊕ Z K √ 2(1 + j) 2 ⊕ Z K 1 + i + j + k 2 .
Since there exist elements λ ∈ Λ satisfying N (λ) = 2 which implies m Λ ( 1 λ ) = 1 2 , we know that M (Λ) ≥ 1 2 . Consequently we have

M (Λ) = 1 2 .
Now, let us come back to our algorithm. A Z-basis of Λ is for instance

1, √ 2, √ 2 2 (1 + i), 1 + i, √ 2 2 (1 + j), 1 + j, 1 2 (1
+ i + j + k), √ 2 2 (1
+ i + j + k)
and, according to this choice of basis, the matrix of the fundamental embedding is

M =              1 √ 2 √ 2/2 1 √ 2/2 1 1/2 √ 2/2 1 - √ 2 - √ 2/2 1 - √ 2/2 1 1/2 - √ 2/2 0 0 √ 2/2 1 0 0 1/2 √ 2/2 0 0 - √ 2/2 1 0 0 1/2 - √ 2/2 0 0 0 0 √ 2/2 1 1/2 √ 2/2 0 0 0 0 - √ 2/2 1 1/2 - √ 2/2 0 0 0 0 0 0 1/2 √ 2/2 0 0 0 0 0 0 1/2 - √ 2/2              .
We use the following cutting-covering of F: if a l = min{x l ; x ∈ F } and b l = max{x l ; x ∈ F} for every 1 ≤ l ≤ 8, we put h l = (b l -a l )/2n l where n 1 = n 2 = 70, n 3 = n 4 = n 5 = n 6 = 30, n 7 = n 8 = 12 and we cover F with parallelotopes P as in [START_REF] Vignéras | Arithmétique des algèbres de quaternions[END_REF] where for every l, C l = a l + (2i l + 1)h l for some 0 ≤ i l ≤ n l -1. Of course, as in the number field case, we only consider parallelotopes that intersect F . Running Algorithm 2 with k = 1/2, X = M • ([-2, 3] ∩ Z) 8 , after four loops of Steps 3, 4, and 5, we finally obtain a convenient digraph composed of nine disjointed simple circuits whose length is 1. This leads to E = {t i = Φ(ξ i ); 1 ≤ i ≤ 9} where 2 j. This leads to M (Λ) = 1/2 as expected.

Of course, and this is the interest of the algorithm, in general we do not have any theoretic argument that allows to establish the value of M (Λ) as in the previous example.

6.2. Some other results. Table 1 gives some of the results that we have obtained. More extensive tables are available from [START_REF] Cerri | Tables of Euclidean minima of orders in totally definite quaternion fields over quadratic fields[END_REF]. In this table, K is a number field, F is a quaternion field and Λ is an order of F . For brevity, the order Λ is described by a quadruple (q 1 , q 2 , q 3 , q 4 ) ∈ F 4 such that Λ = ⊕ 4 i=1 q i Z K , which is always possible in the cases that we have studied. The column Max indicates if Λ is maximal or not (maximal orders were obtained with [START_REF]Handbook of Magma functions[END_REF]), T is the number of critical points in Φ(F ) modulo R = Φ(Λ) and M (Λ) is the Euclidean minimum of Λ. 6.3. Timings. The time required to execute Algorithm 2 obviously depends on the choice of parameters. Indeed, if we choose k > 0 strictly larger than M (Λ), then we will have to start again with a smaller value of k. The timings presented here are for a choice of k smaller than M (K), typically k = M (Λ) -ǫ, where ǫ = 10 -4 .

In Table 2, the orders are the maximal orders corresponding to the quaternion field given in Table 1. It should be noted that no particular efficiency was pursued for Step 6 because it turned out to be short enough in the examples considered.

The times given are total computation times (cpu times). The computations were carried out on an Intel Xeon E5-2680 (2.50GHz) processor with 24 cores.

d F Λ Max T M (Λ) Q( √ 2) -1, -1 K 1, √ 2(1+i) 2 , √ 2

Table 1 .

 1 Euclidean minima of some quaternion fields.

										(1+j) 2	,	yes	9	1/2
										1+i+j+k
											2
	-1, -1 K -3, -√ 2 -2 K		(1, i, j, k) 1, 1+i 2 , j, j+k 2	no yes 36 1	4 41 16
	Q( √ 3)	-1, -1 K	1, 1+i 2 ( √ 3 + 1), 1+j 2 ( √ 3 + 1), 1+i+j+k 2	yes	9	1
		-1, -1 K	1, i,	√	3+j 2 ,	√ 3i-k 2	yes	9	1
	Q( √ 5)	-1, -1 K	1, i, 1+	√ 5+(3+ 4	√	5)i+2j	,	yes	3	1/4
										1+i+j+k
											2
	-1,	√ 5 -7 K	1, i,	√ 5+3 4 (1 + i) + j 2 , √ 5+3 4	+ k 2	yes	6	9/4
	Q( √ 13)	-1, -1 K	1, i, 3+	√ 13 4	+ (5+	√ 13)i 4	yes 24	3/4
								+ j 2 , 1+i+j+k 2
	F		time for loops of			time for	total time
			absorption and units		tests
				(Steps 3 to 5)			(Step 6)
	-1, -1 Q( √ 2) -1, -√ 5 -7 Q( √ 5)			24h16min 7h19min					44min 5 min	25h2min 7h25min
	-1, -1 Q( √ 13)			9d12h36min				1d2h41min	10d15h21min

Table 2 .

 2 Some timings for Algorithm 2.Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux and CNRS, and ANR in accordance to the "programme d'investissements d'Avenir" (see http://www.plafrim.fr).

Acknowledgments

The research of the second author was funded by a "nouveau chercheur" grant by région Auvergne. Parts of this research were conducted during a visit of the second author to Bordeaux, funded by ERC Starting Grant ANTICS 278537. Both authors would like to thank the anonymous reviewer for her/his remarks and corrections.

Algorithm 2: Computation of M (Λ)

Input: Λ, order of quaternion field F over K.

Output: M (Λ) or failure. 1 Cover F with parallelotopes as described above.

2 Choose k > 0, choose a list X of small elements of R: for instance, X ∈ X iff X = M T where T ∈ Z 4n and -B ≤ T i ≤ B + 1 for some bound B ≥ 0 Moreover, the upper bound used in ( 16) is optimal because there exists a vertex S of P which satisfies N R (S -X) = n i=1 A i (P, X).

4.3.

The units test. Let ε = ±1 be the unit of Z K used for this test. Write {P 1 , . . . , P N } for the set of problematic parallelotopes remaining after Step 3, and for every p ∈ {1, . . . , N }, let us denote by D p the center of P p . Let P be one of these parallelotopes, centered in C:

Then ε ⋆ P is a parallelotope centered in C ′′ where C ′′ i+kn = σ i (ε)C i+kn for 1 ≤ i ≤ n and 0 ≤ k ≤ 3 and whose faces are also orthogonal to the canonical axes of R 4n . More precisely

. We first determine a list of elements of R containing all the X ∈ R such that ε ⋆ P -X meets F . Using the same argument as in [START_REF] Cerri | Euclidean minima of totally real number fields: Algorithmic determination[END_REF] we see that such an X is of the form (17) X = X 0 + M (ν 1 , . . . , ν 4n ),


Christine Chataigner

Olivier Barge

Quantitative Approach to the Diffusion of Obsidian in the Ancient Northern Near East

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Obsidian sources located in the northern Near East (Eastern Anatolia, Transcaucasia) were extensively exploited during the Neolithic, Chalcolithic and Early Bronze Age periods, between the 7 th and 3 rd mill. BC. In traditional studies of regional obsidian trade, the direct linear distance from the source is generally considered to be an essential parameter for explaining the quantity of material recovered on the archaeological sites (Renfrew 1977; Peterson / Mitchell / Shackley 1997). However, studies [START_REF] Keller | F. Pernicka / V. Nasedkin, Armenian and Caucasian Obsidian Occurrences and Sources for the Neolithic Trade: Volcanological Setting and Chemical Characteristics[END_REF][START_REF] Blackman | Chemical Characterization of Caucasian Geological Obsidian Sources[END_REF]; Badalyan / Chataigner / Kohl 2004) carried out in Transcaucasia (Fig. 1) show that the factor of direct linear distance often does not apply. Most of the villagers of these prehistoric periods optained obsidian from several sources, and the closest deposit was not necessarily the most desirable. In the whole of the corpus, only 40% of the sites are mainly supplied by the nearest source; 42% exploit the nearest source less than other ones, and 18% do not exploit it at all (Barge / Chataigner 2003). Other factors need to be considered, which is why we have constructed a GIS which takes distances into account, as well as the constraints for movement and the transport of obsidian from the sources to the sites.

Travel Costs

The first stage of our analysis was to evaluate the nature and patterns of travel cost between the Trans caucasian sources of obsidian and the archaeological sites, in order to understand the cost factors actually involved for the prehistoric peoples who sought this material:

Snow is abundant in winter above 2000 m and • blocks access to obsidian sources located at high altitude, but this is only a seasonal constraint, as in summer these highlands become pastures which are frequented by transhumant shepherds.

The vegetation is generally of steppe type and it • was so at the beginning of the Holocene according to palaeo-environmental studies.

The rivers are mainly narrow streams, and many • fords enable crossing; it is the steep relief of the canyons through which the rivers flow, rather than the rivers, which presents difficulties. The topographical element, with elevations often • higher than 3000 m and deep valleys, appears to be the main constraint for travelling. Other, nonenvironmental factors will have also influenced the movement across the ancient landscape: political boundaries, social divisions, cultural taboo or attraction associated with certain sites; however these factors have little if any mark in the archaeological record and therefore can be difficult or impossible to reconstruct (Bell / Wilson / Wickham 2002).

Transport in the 7 th -3 rd mill. BC would have been either on human backs or on the backs of oxen, but not on the backs of equids, which were domesticated locally only in the Early Bronze Age. The effort of the bearer would have been a function of the weight carried. According to ethnographic studies, the pack oxen can carry loads of 50-90 kg, and human bearers loads of 20-80 kg or more. A recent experiment carried out with soldiers provides interesting information on the relations existing between the load of the backpack and the speed of marching, in ensuring that metabolic demands are retained at a bearable level: a) load 50 kg -speed 3.5 km/h; b) load 35 kg -speed 4.5 km/h; c) load 20 kg -speed 5.5 km/h [START_REF] Scott | [END_REF].

From this experiment, we can infer for a reasonably trained walker, an average speed on the flat of 5 km/h, for a load between 25-30 kg and a walking time of 7 to 8 hours per day. Even if it is difficult to extrapolate from the objects found in an excavation the total weight of obsidian which could be brought to the site, such a quantity of obsidian (30 kg) is already a considerable amount; for example, the neolithic site of Aratashen, in the Ararat plain (Armenia), has yielded more than 18,000 obsidian artefacts for a total weight of about 90 kg [START_REF] Badalyan | New data on the late prehistory of the Southern Caucasus. The excavations at Aratashen (Armenia): preliminary report[END_REF].

The main travel cost being the topography, the cost surface of the GIS has been defined by the time needed for a walker depending on the slope and the distance. Two formulas have been proposed, with results comparable for the slopes less than 25° (near the maximal gradient accepted by mountain trekkers):

(1) t = d / 6e -3.5 |s + 0.05| (Gorenflo / Gale 1990), where t = time in hours, d = distance in km and s = slope, calculated as vertical change divided by horizontal change; e is the exponential function (2) t = (0.031 s 2 -0.025 s + 1) * d/5 (Eastman 1999;Schneider / Robbins 1996), where t = time in hours, d = distance in km and s = slope angle in degrees. These formulae are anisotropic, since the movements uphill and downhill will produce different results. However, in the particular case of obsidian procurement, the deposits are located at high altitude, with an ascending outward journey, but with very little load, while the return journey benefits in the descent, but with a much heavier load (obsidian blocks) and fatigue. Also in the hypothesis in which pack animals participated in the transport of the obsidian, the speed of descent would have been limited by the load carried and the difficulties of the terrain. Therefore, wishing to take into account the equivalence of effort which can exist between an ascent with no load and a descent with a load, we used Eastman's formula with an unsigned slope angle, in order to calculate full cost uphill, no cost crosshill and full cost downhill.

The calculations were made for speeds corresponding to a loaded person walking (5 km/h), the results being not much different for draft animals. The ethnographic data indicate that the number of walking hours per day, which would be about 8 hours for a human bearer, is reduced to about 5 hours for draft animals, as they need to stop to graze.

The costweighted distance analysis functions enabled calculation on the one hand of the times accumulated to go from a given obsidian source to the archaeological sites and on the other hand to visualize the best route (or leastcost route) to take. The diffusion from a source of obsidian can be ex- pressed by isolines, which enable an understanding of the actual time necessary to travel from the source to the sites, and to compare for a given region the diffusion of obsidian from different sources (Fig. 2).

Overland Transport: Transcaucasia

In the case of Transcaucasia, important variability can be noticed in the time calculated between sources and sites. In order to better understand where this variability comes from, let us examine a source, Arteni, on a finer scale (isochrones of 2 hours) (Fig. 3a). We notice that the frequency of the sites which have received obsidian from this source as well as the percentage of this obsidian in their material decreases dramatically above the isochrone of 15 hours.

The diffusion of the obsidian of Arteni can be represented by a histogram, which shows, by time ranges of 7 hours (or one day walk), the percentage of sites present in this range and the average per-centage of obsidian which they use (Fig. 3b). The threshold at 14/15 h is clearly visible.

Arteni is situated at the periphery of the Ararat plain, at an altitude of about 1500 m, and is ge nerally accessible all year round. Two other important obsidian sources, Gutansar and Hatis, are similarly situated (Badalyan / Chataigner / Kohl 2004, 450-453). The histograms of these three sources show a great similarity in the manner of diffusion (Fig. 3b,c). We notice at first a threshold at 15 hours (about 2 days walk). This threshold could correspond to the maximum acceptable time of direct procurement for the villagers.

For the villages which are situated between 15 to 42 hours from the source, the obsidian is present in reduced quantity: this probably signifies redistribution. Examination of the leastcost routes between the source of Arteni and these farther sites shows that they pass through intermediate sites, from which obsidian could have been passed on.

Other groups of obsidian sources (Geghasar, Sevkar) (Badalyan / Chataigner / Kohl 2004, 453- 455) are found on high plateaus at altitudes of 2500-3000 m. The histograms based on the time needed to access the sources indicate again the existence of a threshold at about 14/15 hours, a threshold which concerns here the quantity of obsidian and not that of the sites, which remain rare (Fig. 4a). The map of the modelled pathways shows that this first group of sites is located on the west shore of Lake Sevan, at an elevation of ca. 2000 m, under severe climatic conditions (Fig. 4b). These villages obtained their obsidian almost exclusively from the Geghasar source. At 15 to 42 hours from the source, we notice a second group of sites, which, in the case of Geghasar, obtained up to 50% of obsidian from this source (Fig. 4a). If we examine the map of routes, we notice that, contrary to what has been observed for Arteni, these distant sites are linked directly to the source (Fig. 4b). It is probable that these villages, established in the valley of the Arax, were those of semi transhumants, who drove their herds in summer to the high plateaus, like nowadays.

The ethnographic sources provide evidence of the practice of this custom in past centuries and the absolute necessity to act in this way, because of the very arid and very hot summer climatic conditions in the valley (Mkrtumyan 1974).

A third group of sites is found even farther away, up to 70 hours by foot from the source (about 10 days for a human bearer and 14 days for pastoralists with their herd), and we notice in the case of Sevkar that the percentage of obsidian coming from the source can be high (Fig. 4a). Here again study of the routes shows that some of these sites are directly linked to the source (Fig. 4b). This suggests the hypothesis of "longdistance" transhumants, whose herds grazed yearround in pastures, in the winter in the mild climate steppes of Azerbaijan and in the summer on the high plateaus overlooking Lake Sevan. These transhumants brought down a noticeable quantity of obsidian, probably thanks to pack animals.

The example of Geghasar shows us also that several modes of distribution can exist for the same source, as there also seems to be evidence of redistribution from near settlement to near settlement for sites established along the Agstev river and in the Kura valley (Fig. 4b).

For the Sevkar source, we find three similar patterns of direct acquisition: 1) by villagers in the vicinity, in the Vorotan valley; 2) by semitranshumants, settled along the shore of the Lake Sevan or in the Araxes valley; 3) by longdistance transhumants, coming from the steppes of Azerbaidjan or from the lake Urmiah basin. This hypothesis of a significant transhumance suggested by the GIS is confirmed archaeologically by the excavation of the site of Godedzor, situated in the southeastern part of the Lesser Caucasus (Chataigner et al. in press).

This Late Chalcolithic settlement is situated at an altitude of about 1800 m in a region covered by snow from November to April. This settlement has yielded an enormous quantity of material, but very few architectural remains, suggesting that this village was deserted during the snowy winter. Almost all of the lithic industry of Godedzor is in obsid- ian. This material was obtained from the deposits of Sevkar and Bazenk, situated on the high plateaus which dominate the region. The inhabitants of Godedzor also possessed painted pottery which originated in northwestern Iran, specifically in the basin of Lake Urmiah. Moreover, the herd was composed solely of sheep and cattle and several phalanges and vertebras of oxen presented deformations due to carrying heavy loads. All these elements suggest that Godedzor was a seasonal village and that its inhabitants were pastoralists who came from the basin of Lake Urmiah, to which they returned in the autumn with obsidian blocks from the Sevkar and Bazenk deposits.

Water Transport: Lake Van

In the mountains of the northern Near East and Trancaucasia, the rivers are not navigable. On the other hand the three large lakes of the region could have been crossed by boat.

There is evidence for navigation in the ancient Near East at least since the 9 th millennium, as the site of Shillourokambos on the island of Cyprus received obsidian from sources in Cappadocia (Briois / Gratuze / Guilaine 1997). There are no remains of Neolithic boats in the Near East, but it is probable that they were carved in tree trunks, like the one found on a Neolithic site of the 6 th millennium in Italy (Lake Bracciano). To test the viability of this logboat, a Czech team built a faithful copy of the archaeological example, and made a voyage along the Mediterranean seacoast, with provisions and 100 kg of obsidian from Lipari (whose blocks were found on the Italian site) (Tichy 2000).

The replica built according to the find from Bracciano proved to be very good on the sea. The transport of obsidian was optimal, since the load was used as ballast. The average speed on the 800 km of the route was 4 km/h which meant mostly covering a distance of more than 50 km in a day.

The Nemrut Dag volcano, whose obsidian was widely diffused throughout the Near East from the beginning of the Neolithic, dominates the western extremity of Lake Van (Fig. 5). Its obsidian was diffused in particular towards the east, to Tilki Tepe; thus it seemed interesting to test whether transport by boat on Lake Van could have played a role in this pattern. We consequently carried out two calculations with the costweighted distance function, the first considering that the lake could have been crossed at the speed of 4 km/h, and the second considering that the lake had to be circumvented at the speed of man walking (5 km/h).

In the first case, the crossing would have taken about 25 hours. The accounts of travellers in past centuries refer to the crossing by boat from Tatvan to Van as not very dangerous and quite rapid (about 24 hours) (Yerasimos 1991), which corresponds to an average speed of 4.2 km/h, which agrees with the logboat experiment.

If the lake had to be skirted, the route along the north shore would have taken 45 hours. Reports of travellers of earlier centuries confirm the length of this route, which required 8 days by horse, and stress that the southern route is geographically shorter but more difficult. According to our assumptions, the total time by the southern route is in fact equivalent to that for the northern route. The quantity of obsidian found at Tilki Tepe is considerable, so the hypothesis of transport by boat, being more rapid and less tiring, is certainly plausible.

Conclusion

In this initial stage of the study of the circulation of obsidian in the ancient Near East using a GIS, we have especially sought to verify the validity of the modelling. The systematic comparison between the hypotheses suggested by the GIS and other relevant data (archaeological, ethnographical or experimental data etc.) clarifies and supports the results of the modelling achieved using the GIS.

Fig. 1 .

 1 Fig. 1. Location of the obsidian sources in the northern Near East (Transcaucasia, Eastern Anatolia).

Fig. 2 .

 2 Fig. 2. Variability of the diffusion of the Transcaucasian obsidian sources: the examples of (a) Arteni, (b) Geghasar and (c) Chikiani.

Fig. 3 .

 3 Fig. 3. Diffusion of sources of obsidian situated at the periphery of the plain of Ararat: (a) Arteni: quantitative distribution of obsidian. (b)-(c) Arteni, Hatis and Gutansar: variations in the quantity of obsidian and the number of sites in relation to the time needed for access.

Fig. 4 .

 4 Fig. 4. Diffusion of the sources of obsidian situated at more than 2500 m in altitude: (a) Geghasar and Sevkar: histograms of the diffusion of obsidian in relation to the time needed for access. (b) Geghasar: leastcost pathways between the source and the sites.