Vincent Despré
email: vincent.despre@loria.fr

Loïc Dubois
email: loic.dubois@ens-lyon.fr

Benedikt Kolbe
email: benedikt.kolbe@physik.hu-berlin.de

Monique Teillaud
email: monique.teillaud@inria.fr

Experimental analysis of Delaunay flip algorithms on genus two hyperbolic surfaces

Guided by insights on the mapping class group of a surface, we give experimental evidence that the upper bound recently proven on the diameter of the flip graph of a closed oriented hyperbolic surface by Despré, Schlenker, and Teillaud (SoCG'20) is largely overestimated.

To this aim, we develop an experimental framework for the storage of triangulations. We show that the computations with algebraic numbers can be overcome by proving a density result on rationally described genus two surfaces, and we propose ways to generate surfaces that are meaningful for the experiments.

Introduction

It was recently proven that the geometric flip graph of a closed oriented hyperbolic surface is connected [START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF]. A Delaunay flip algorithm can thus transform any input geometric triangulation T , i.e., a triangulation whose edges are embedded as geodesic segments only intersecting at common endpoints, into a Delaunay triangulation. This is particularly useful in practice as a crucial preprocessing step to computing Delaunay triangulations on a surface: it transforms a "bad" representation of a surface, e.g., by a very elongated fundamental domain, to a "nice" representation by a Delaunay triangulation with only one vertex.

An upper-bound on the number of flips was proven [START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF]Theorem 19]: C h • ∆(T ) 6g-4 • n 2 , where C h is a constant, ∆(T ) is the diameter of T , g is the genus of the surface, and n is the number of vertices. The diameter ∆(T ) is the smallest diameter of a fundamental domain that is the union of lifts of the triangles of T in H (note that this is not the diameter of the surface, which is independent of the representation). Computing it algorithmically seems difficult, however for a triangulation with only one vertex (thus with 4g -2 triangles) some bounds are easily derived:

L T ≤ ∆(T ) ≤ ∆(F ) ≤ L T .(4g -2)
, where L T denotes the maximal length of an edge of T and F ⊂ H is any fundamental domain made of lifts of the triangles of T . From these bounds we see that ∆(T ) cannot differ too much from the diameter of any such F : in the case of a genus two surface they only differ by a factor of at most six. In the experiments, we will thus use the domain that naturally appears.

In this paper, we experimentally study the dependence of the number of flips on ∆(T ) (Section 7), for surfaces of genus two. We suspect that the factor ∆(T ) 6g-4 is largely overestimated.

It is derived from the number of paths of bounded length on a surface. Intuitively, for a length bounded by L, it roughly amounts to the volume of the ball of diameter L, so, it is exponential in L; if only simple paths are considered, this number reduces to L 6g-4 [START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF], but there is no reason why the flip algorithm would use all the simple paths shorter than L instead of going straight.

More formally, our expectation on the dependence in ∆(T ) is based on insights on the structure of the mapping class group (Section 2.3).

To perform experiments, we set up a framework consisting of various tools. In Section 3, we present a data structure for triangulations of surfaces, which supports flips; it relies on the representation of genus two surfaces by octagons in H (Section 2.4). Not surprisingly, arithmetic issues quickly arise, as algebraic numbers are involved in the description of the octagons (Section 4.1). We overcome them by proving a density result on rationally described octagons (Section 4.2), which allows us to restrict to rational numbers in our experiments.

The generation of input surfaces and triangulations is far from trivial; it is a non-negligible part of our work (Section 5). We obtain surfaces with a large diameter by twisting the abovementioned octagons (Section 5.1).

In Section 6, we run experiments comparing strategies on the sequence of edge flips, and conclude that the naive strategy is close to being the best one in practice. We adhere to it for our main experiments that study the dependence of the number of flips on ∆(T ).

The way we conduct these experiments in Section 7 is inspired by previous work by Mark

Bell [START_REF] Mark | Simplifying triangulations[END_REF] who studied flips in a topological setting. We focus on triangulations having only one vertex, both because the dependence on the number of vertices is easily seen, and because inserting a lot of points would rather be done by Bowyer's incremental algorithm [START_REF] Iordanov | Implementing Delaunay triangulations of the Bolza surface[END_REF][START_REF] Despré | Representing infinite hyperbolic periodic Delaunay triangulations using finitely many Dirichlet domains[END_REF], inspired from previous work in the flat case [17]. Quite surprisingly, in practice, we observe a behavior that is only expected asymptotically.

Background

Hyperbolic surfaces

All the surfaces considered in this paper are closed (connected, compact and without boundary), oriented, and hyperbolic. Consider such a hyperbolic surface S of genus g ≥ 2 and the underlying topological surface S g . Given a hyperbolic structure h on S, associated to a metric of constant curvature -1, the surface S = (S g , h) is isometric to the quotient H/G, where H is the hyperbolic plane and G is a (non-Abelian) discrete subgroup of the isometry group of H isomorphic to the fundamental group π 1 (S g ).

The universal cover of S is isometric to H equipped with a projection ρ : H → S that is a local isometry. The group G acts on H, so that for any p ∈ S, ρ -1 (p) is an orbit under the action of G. A lift p of a point p ∈ S is one of the elements of the orbit ρ -1 (p).

We use the

Triangulations and flips on hyperbolic surfaces

A topological triangulation of a hyperbolic surface S is any embedding of an undirected graph with a finite number of vertices onto S such that each resulting face is homeomorphic to an open disk and is bounded by exactly three distinct edge-embeddings. Observe that this graph may have loops or multiple edges, and recall that the terms embedding and embedded subsume that edges only intersect at common vertices. A geometric triangulation is a topological triangulation of S whose edges are embedded as geodesic segments [START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF]. All triangulations considered in this paper are geometric, so we just use the term triangulation. For any triangulation T of S, the lift

T

Mapping class group

We use the same notation as Maher [START_REF] Maher | Random walks on the mapping class group[END_REF] and refer to his paper for details.

The set Mod(S g ) of all homeomorphims (up to isotopy) of a topological surface S g is called the mapping class group of S g . Following Thurston's classification [START_REF] Fathi | Thurston's Work on Surfaces (MN-48)[END_REF], Mod(S g ) contains three kinds of elements: the periodic homeomorphims, which are of finite order and are not useful for our purposes; the reducible ones, which fix at least one curve on S g ; and the so-called pseudo-Anosov homeomorphims, also known as the hyperbolic elements of Mod(S g ).

Dehn twists (Figure 1, Left) are typical reducible elements, as they fix all the curves that do not intersect the curve used for twisting. A Dehn twist by a curve c at most adds to the length of a curve a constant that depends on the number of times the curve intersects c. A pseudo-Anosov element at most multiplies the length of a curve by a constant factor. Mod(S g ) can be generated by a finite set of Dehn twists [START_REF] Farb | A primer on mapping class groups, volume 49 of Princeton Mathematical Series[END_REF]. The composition of generators or their inverses in a random order can be interpreted as a random walk in Mod(S g ): such a walk reaches pseudo-Anosov elements with asymptotic probability 1 [START_REF] Maher | Random walks on the mapping class group[END_REF]. However, this asymptotic result does not a priori describe the local structure of Mod(S g ).

Admissible symmetric octagons

The Teichmüller space T M 2 of the topological surface S 2 is the set of all the hyperbolic structures (up to isotopy) that can be associated to S 2 . It admits various parametrizations. The most commonly used, though not well adapted to our needs, is the set of Fenchel-Nielsen coordinates [START_REF] Hubbard | Teichmüller Theory and Applications to Geometry, Topology, and Dynamics[END_REF]Section 7.6]. Here, we use a less usual set of parameters introduced by Aigon-Dupuy, Buser et al. [START_REF] Aigon-Dupuy | Hyperbolic octagons and Teichmüller space in genus 2[END_REF], who proved that any surface of genus 2 has a fundamental domain that is an octagon in D. This versatile representation allows us to easily construct and manipulate such surfaces in our experiments. In this section we recall some definitions and results of the original paper [START_REF] Aigon-Dupuy | Hyperbolic octagons and Teichmüller space in genus 2[END_REF], following its notation. Given j ≥ 3 and complex numbers z 1 , . . . , z j ∈ D in convex position, G[z 1 , . . . , z j ] denotes the hyperbolic polygon whose vertices are z 1 , . . . , z j in this order. Given a compact subset X ⊂ D, A(X) is the hyperbolic area of X. Given z ∈ C, we denote by Re [z] and Im [z] its real and imaginary parts, respectively, by z its conjugate, and by |z| its modulus; i denotes a root of -1.

Let arg z ∈ [0, 2π[ denote the argument of a point z = 0 C . Given z 0 , z 1 , z 2 , z 3 ∈ D \ {0 C }, the 4-tuple (z 0 , z 1 , z 2 , z 3 ) is valid if 0 = arg z 0 < arg z 1 < arg z 2 < arg z 3 < π; the hyperbolic octagon P[z 0 , z 1 , z 2 , z 3 ] is then defined as G[z 0 , z 1 , z 2 , z 3 , -z 0 , -z 1 , -z 2 , -z 3 ]
. Such a hyperbolic octagon is called a symmetric octagon. The interior angles of a symmetric octagon cannot be greater than π. If moreover A(P[z 0 , z 1 , z 2 , z 3 ]) = 4π, then P[z 0 , z 1 , z 2 , z 3 ] and the 4-tuple (z 0 , z 1 , z 2 , z 3 ) are called admissible. Each surface of genus 2 can be obtained by identifying the opposite sides of an admissible symmetric octagon [START_REF] Aigon-Dupuy | Hyperbolic octagons and Teichmüller space in genus 2[END_REF]. Observe that the eight vertices of the octagon correspond to the same point on the surface.

A valid 4-tuple (z 0 , z 1 , z 2 , z 3 ) is admissible if and only if Im

3 k=0 (1 -z k z k+1 ) = 0 [1, Lemma 3.2]
. The authors establish this condition after proving a preliminary result that we will reuse: for any two points z, z

∈ D \ {0 C } if 0 ≤ arg z ≤ arg z ≤ π then [1, Appendix (A7)] 2 arg(1 -zz ) = A(G[0 C , z, z ]). (1)
An admissible 4-tuple can be constructed as follows [START_REF] Aigon-Dupuy | Hyperbolic octagons and Teichmüller space in genus 2[END_REF]Section 3].

Start with z 1 , z 2 , z 3 ∈ D satisfying 0 < arg(z 1 ) < arg(z 2 ) < arg(z 3 ) < π. Abbreviate u = (1 -z 1 z 2 )(1 -z 2 z 3 ), a = Im [-uz 1 z 3 ], b = Im [u(z 3 -z 1 )], and c = Im [u]. Assume a + b + c < 0 and let z 0 = 2c -b+ √ b 2 -4ac .
Then (z 0 , z 1 , z 2 , z 3 ) is an admissible 4-tuple.

From now on indices are modulo 8. Let us consider an admissible 4-tuple (z 0 , z 1 , z 2 , z 3 )

and define z l+4 = -z l for every l ∈ {0, 1, 2, 3}. For k ∈ {0, . . . , 7}, there exists a unique orientation preserving isometry τ k : D → D satisfying τ k (z k+5 ) = z k and τ k (z k+4 ) = z k+1 :

the isometry τ k maps a side of P[z 0 , z 1 , z 2 , z 3 ] to the opposite side of P[z 0 , z 1 , z 2 , z 3 ]. Define

ω k = z k (1 -|z k+1 | 2 ) + z k+1 (1 -|z k | 2 ) 1 -|z k z k+1 | 2
and note that |ω k | < 1; the isometry τ k is then given by

τ k (z) = (z + ω k )/(ω k + 1) for every z ∈ D [1, Lemma 4.1]. Observe that τ k+4 = τ -1 k .

Representation of triangulations

In this section we describe our data structure for representing triangulations (Section 3.1) and we sketch how it is maintained through flips (Section 3.2).

Data structure

Although an ad hoc data structure was previously proposed for flipping triangulations [START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF], we choose to use combinatorial maps, which are commonly used to represent graphs embedded on a surface. We refer the reader to the literature for a formal definition [START_REF] Mohar | Graphs on Surfaces[END_REF]Section 3.3]. The data structure we use offers a representation of the triangulation that intrinsically lies on the surface, while the earlier data structure [7, Section 4.1] stuck to specific representatives of all vertices and faces of the lifted triangulation in the universal cover. See Appendices A and B.1 for details on this section.

For our experiments, we use the flexible implementation of combinatorial maps that is publicly available in cgal [START_REF] Damiand | CGAL User and Reference Manual[END_REF]. The dart, also known as flag, is the central object in a combinatorial map: it gives access to all incidence relations of an edge of the graph (Figure 2). In our setting a combinatorial map can be thought of as a halfedge data structure. The geometric information for the triangulation is stored by adding a cross-ratio for each edge.

Recall that the cross-ratio of four pairwise-distinct points in H represented by Note that in our experiments, the lifts in D are only used to calculate the cross-ratios of a given input triangulation T ; they are ignored during the flips, thus preserving the property that the data structure only considers the embedding of the triangulation on the surface. However, in order to be able to recover a lift in D in the end, e.g., for drawing a representation in D of the final Delaunay triangulation, we need to maintain an anchor during flips. The anchor A = (δ, a 1 , a 2 , a 3 ) consists in a dart δ, chosen arbitrarily, together with a triple (a 1 , a 2 , a 3 ) of the vertices of a lift of the face containing δ.

z 1 , z 2 , z 3 , z 4 ∈ D is the complex number [z 1 , z 2 , z 3 , z 4 ] = (z 4 -z 2 )(z 3 -z 1 ) (z 4 -z 1 )(z 3 -z 2 ) [ 3
A triangulation T is thus represented by (M, F, A), where M is the combinatorial map, F the map that associates a cross-ratio to each edge of M , and A = (δ, a 1 , a 2 , a 3 ) is the anchor.

Flipping an edge

In this section, we quickly sketch how the data structure is maintained through an edge flip.

First we modify the combinatorial map, then we update the anchor, and we finally update the cross-ratios. Some details and the pseudo-code are given in Appendix B.2.

Performing a flip in the combinatorial map is a straightforward use of the functionalities given by the cgal package [START_REF] Damiand | CGAL User and Reference Manual[END_REF]. The triangulation obtained from T after flipping an edge e is denoted by T . By definition, the dart δ of the anchor A belongs to the face t of T represented by a lift t = (a 1 , a 2 , a 3 ) in D. If t does not contain e then A is not modified by the flip. However, if e is an edge of t then t will not belong to T and we must update A. A lift e of e incident to t is replaced by e when e is flipped. The anchor is updated so that it represents one of the two faces incident to e in T .

Finally, the cross-ratios must be retrieved. Only the cross-ratios of the at most 5 edges of the two triangles forming the quadrilateral whose diagonal is to be flipped must be updated.

Their values after the flip are expressed in terms of their values before the flip (see Lemma 5 in Appendix A).

Solving arithmetic issues

The construction recalled in Section 2.4 shows that the real and imaginary parts of the complex numbers involved when defining surfaces are in general algebraic numbers. Efficiency issues when computing with algebraic numbers have been known for decades. More recently, they appeared when constructing Delaunay triangulations of hyperbolic surfaces [START_REF] Iordanov | Implementing Delaunay triangulations of the Bolza surface[END_REF][START_REF] Ebbens | Delaunay triangulations of generalized Bolza surfaces[END_REF], showing that the hope to get effective software was restricted to very simple cases. In Section 4.1 we describe a simple experiment on the Bolza surface illustrating that these arithmetic issues are actually prohibitive in practice for the Delaunay flip algorithm in the sense that they imply unreasonable running times.

We show in Section 4.2 that any surface of genus 2 can be approximated by a surface described by rational numbers. It is straightforward to check that the computations made during a Delaunay flip algorithm only use the four basic operations +, -, •, / (see Section 3.2 and Appendix B.2). Thus if the input surface is represented by rational numbers all numbers arising throughout the algorithm stay rational. This fact allows us to run extensive experiments.

Issues when using algebraic numbers

Let c k = exp iπ 2k-1 8 2 1/4
, k ∈ {0, . . . , 7} be the vertices of a regular hyperbolic octagon in D;

identifying the opposite sides of this octagon gives a surface of genus 2 known as the Bolza surface. Consider the triangulation T 0 of the octagon shown in Figure 3. Identifying in T 0 the edges corresponding to opposite sides of the octagon yields a triangulation T of the Bolza surface. Let e 0 , . . . , e 4 be the edges of T corresponding to the edges e 0 , . . . , e 4 of T 0 . The algebraic

numbers cos π 8 = 2 + √ 2 2
and sin

π 8 = 2 - √ 2 2
naturally appear when computing the cross-ratios R T (e l ) = [c 0 , c l+1 , c l+2 , c l+3 ], l ∈ {0, . . . , 4}.

Density of the rationally described surfaces

is the hyperbolic distance in D. Definition 1. We say that a 4-tuple (z 0 , z 1 , z 2 , z 3 ) is rational if z k ∈ Q + iQ for every k ∈ {0, 1, 2, 3}. A rationally described surface is a surface obtained from a rational admissible 4- tuple (z 0 , z 1 , z 2 , z 3 ) by identifying the opposite sides of P[z 0 , z 1 , z 2 , z 3 ].
Theorem 2. Let (z 0 , z 1 , z 2 , z 3 ) be an admissible 4-tuple and ε > 0. There exists a rational admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ) such that ∀k ∈ {0, 1, 2, 3}, z k ∈ B(z k , ε).

(G[-z 0 , z 0 , z 1 , z 2 , z 3 ]) = 2π
(equivalent to the condition A(P[z 0 , z 1 , z 2 , z 3 ]) = 4π given in Section 2.4). Consider Figure 4.

If ε is small enough, then (z 0 , z 1 , z 2 , z 3 ) is valid. We will now show that if each z k is "close enough" to z k for every k, we can replace z 3 by a point U in B(z 3 , ε) ∩ (Q + iQ) so that the area condition is satisfied. More details on the construction can be found in Appendix C.2.

To do so we first define an isometry f :

D → D in the Poincaré disk: f (z) = z + z 0 z 0 z + 1 . Observe that f (-z 0 ) = 0 C . Since f and f -1 both map D ∩ (Q + iQ) to some subset of D ∩ (Q + iQ) our problem reduces to replacing f (z 3 ) by an element V of B(f (z 3 ), ε) ∩ (Q + iQ) satisfying A(G[0 C , f (z 0 ), f (z 1 ), f (z 2 ), V ]) = 2π. Indeed, by setting U = f -1 (V ) we obtain U ∈ B(z 3 , ε) ∩ (Q + iQ) and A(G[-z 0 , z 0 , z 1 , z 2 , U ]) = A(G[0 C , f (z 0 ), f (z 1 ), f (z 2 ), V ]) = 2π.
To find such a point V , we define a polynomial P ∈ Q[X] by setting

P (X) = Im (1 -f (z 0 )f (z 1 ))(1 -f (z 1 )f (z 2 ))(1 -Xf (z 2 )f (z 3 )) .
Observe that the degree of P is at most 1, and thus P (X) = (P (1) -P (0))X + P (0). We first show that if we choose z k close to z k for every k ∈ {0, 1, 2, 3}, P (1) is close to 0 and P (0) close to κ > 0. Since 0 = arg f (z 0 ) < arg f (z 1 ) < arg f (z 2 ) < arg f (z 3 ) < π we can apply Equality (1)

0 C = f(-z 0 ) z 0 = f(0 C ) -z 0 z 3 z 2 z 1 f(z 3 ) f(z 2 ) f(z 1 ) f(z 0 )
Figure 4: Illustration of the proof of Theorem 2 .

and obtain

arg 1 -f (z 0 )f (z 1 ) 1 -f (z 1 )f (z 2 ) 1 -f (z 2 )f (z 3 ) = arg 1 -f (z 0 )f (z 1 ) + arg 1 -f (z 1 )f (z 2 ) + arg 1 -f (z 2 )f (z 3 ) = 1 2 [A (G[0 C , f (z 0 ), f (z 1 )]) + A (G[0 C , f (z 1 ), f (z 2 )]) + A (G[0 C , f (z 2 ), f (z 3 )])] = 1 2 A (G[0 C , f (z 0 ), f (z 1 ), f (z 2 ), f (z 3 )]) = 1 2 A (G[-z 0 , z 0 , z 1 , z 2 , z 3 ]) .
By observing that every expression in between the equalities belongs to [0, 2π[ we see that those equalities are indeed equalities and not only congruences modulo 2π. By choosing z k close to z k for every k ∈ {0, 1, 2, 3} we make the last expression approach

1 2 A(G[-z 0 , z 0 , z 1 , z 2 , z 3 ]) = π,
which makes P (1) tend to 0. Similarly, we obtain

arg 1 -f (z 0 )f (z 1 ) 1 -f (z 1 )f (z 2 ) = 1 2 A (G[-z 0 , z 0 , z 1 , z 2 ]) .
By choosing z k closer and closer to z k , the last expression tends to 1 2 A(G[-z 0 , z 0 , z 1 , z 2 ]) which is not congurent to 0 modulo π. Thus P (0) is close to some constant κ > 0, whence we can assume that P (1) = P (0).

To construct V set λ = P (0) P (0)-P (1) and let V = λf (z 3 ); we have both V ∈ Q + iQ and P (λ) = 0. We proved that P (1) tends to 0 and that P (0) tends to κ > 0 so λ tends to 1 and V tends to f (z 3 ). Finally, observe that

P (λ) = 0 implies A(G[0 C , f (z 0 ), f (z 1 ), f (z 2 ), V ]) = 2π by
Equality (1).

Remark 3. This theorem implies the density of the hyperbolic structures corresponding to rational admissible 4-tuples in T M 2 with its canonical topology. However, a proof would go beyond the scope of this paper and would be quite technical.

The generation of the input for the experiments

We generate input for the Delaunay flips algorithms by triangulating admissible loosely-symmetric octagons; the latter are defined in Section 5.1 as a slight extension of the admissible symmetric octagons (Section 2.4). We generate such octagons whose vertices are represented in D by complex numbers with rational real and imaginary parts. To conduct experiments, we need to generate a large number of triangulations with a large diameter. An effective approach consists in generating octagons with a small diameter and twisting them (Section 5.1) many times to obtain octagons with a very large diameter. This approach has two advantages. Firstly it can generate triangulations with different diameters but lying on a single surface, hence eliminating any dependency on the choice of the surface if needed. Secondly thanks to this approach we will also study the dependency of the number of flips on those twists.

The twists of admissible loosely-symmetric octagons

We say that a hyperbolic octagon P is loosely-symmetric if the opposite sides of P are isometric and the opposite interior angles of P are equal. If moreover A(P ) = 4π then we call P admissible.

Clearly Let G[z 0 , . . . , z 7 ] be an admissible loosely-symmetric octagon. We will consider the Dehn twists along the axes of its side-pairings (Figure 1). These twists generate a subgroup of Mod(S 2 ) (Section 2.3), which contains non-reducible elements of Mod(S 2 ) since the generators do not all fix a common curve. Thus, this subgroup contains pseudo-Anosov elements [START_REF] Ivanov | Subgroups of Teichmüller modular groups[END_REF].

For every k ∈ {0, . . . , 7} we denote by τ k the orientation preserving isometry of D that satisfies τ k (z k+5 ) = z k and τ k (z k+4 ) = z k+1 . We fix t ∈ {0, . . . , 7}. For k ∈ {0, . . . , 7} let

z k = τ t (z k ) if k -t ∈ {1, 2, 3, 4} mod 8, z k otherwise.
By the Gauss-Bonnet formula, the interior angles of G[z 0 , . . . , z 7 ] sum up to 2π; since opposite interior angles are equal, each interior angle is at most π. Thus the geodesic segment between z t+1 and z t+5 is contained in G[z 0 , . . . , z 7 ] and cuts the polygon into the two interior disjoint pentagons

P 1 = G[z t+1 , z t+2 , z t+3 , z t+4 , z t+5 ] and P 2 = G[z t+5 , z t+6 , z t+7 , z t , z t+1
]; the intersection of P 1 and P 2 is the segment between z t+5 and z t+1 and the two pentagons are isometric. Similarly, τ t (P 1 )

and P 2 are interior disjoint, they intersect on the segment between z t = z t and z t+4 = z t+1

and their union is G[z 0 , . . . , z 7 ]. It follows that G[z 0 , . . . , z 7 ] is an admissible loosely-symmetric octagon; the surface that it defines is isometric to the one defined by G[z 0 , . . . , z 7 ] as both can be obtained by the same identification of the sides of P 1 and P 2 (P 1 and τ t (P 1 ) being isometric).

We say that (z 0 , . . . , z Let us denote by (τ k ) 0≤k≤7 the isometries defined for z 0 , . . . , z 7 , in the same way as (τ k ) 0≤k≤7 above. By definition of the t-twist, the following holds for every k ∈ {0, . . . , 7}

τ k =    τ t • τ k if k -t ∈ {1, 2, 3} mod 8, τ k • τ -1 t if t -k ∈ {1, 2, 3} mod 8, τ k if k = t mod 4.
For a word t = t 1 • • • t m , we define the t-twist as the composition of the t k -twists, k = 1, . . . , m, in this order. We pick t 1 , . . . , t m in {0, . . . , 3} m instead of {0, . . . , 7} m to only consider the generators without their inverses and obtain large diameters as quickly as possible.

The generation of triangulations

The input surfaces and triangulations are constructed in four steps. We refer to Appendix D for details.

[step 1] We construct an initial rational admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ).

[step 2] We choose n p ≥ 0 and construct points (p 1 , . . . , p np ) ∈ (Q + iQ) n lying within the closure of P[z 0 , z 1 , z 2 , z 3 ].

[step 3] We choose m ≥ 0 and a sequence t 1 • • • t m of twists.

[step 4] From the 4-tuple (z 0 , z 1 , z 2 , z 3 ), the points (p 1 , . . . , p np ), and the sequence

t = t 1 • • • t m ,
we construct a representation (M, F, A) of an input triangulation T . This is (roughly) done by t-twisting P[z 0 , z 1 , z 2 , z 3 ], triangulating the octagon resulting from the twists, and inserting the points in the faces of the resulting triangulation. Together with the point corresponding to the vertices of the octagon, the triangulation T has n = n p + 1 vertices.

Sections 6 and 7 refer to these four steps.

Step 1 was applied a thousand times to construct the 1,000 rational admissible 4-tuples Q 1 , . . . , Q 1,000 ; in our experiments we consider the first n q 4-tuples. In step 3 we constructed 10,000 random sequences of twists denoted by S 1 , . . . , S 10,000 , each of length 10, of which some of the experiments use the first n s sequences. The values of n p , n q , and n s are specified in the description of each experiment.

Technicalities for steps 1, 2, and 4 are deferred to Appendix D. We only elaborate on step 3 here. Consider a sequence of m twists represented by the word t = t 1 . . . t m (see Section 5.1).

We will study two kinds of sequences in the experiments of Sections 6 and 7:

• A power sequence is represented by a word u m for some u ∈ {0, . . . , 3}.

• In a random sequence, t 1 , . . . , t m are chosen uniformly and independently in {0, . . . , 3}.

It appears in practice that the length of a random sequence has a stronger impact on the computations than the length of a power sequence. When twisting, we update an 8-tuple (z 0 , . . . , z 7 ) ∈ (Q+iQ) 8 corresponding to an admissible loosely-symmetric octagon G[z 0 , . . . , z 7 ],

together with the orientation preserving isometries (τ k ) 0≤k≤7 identifying its opposite sides. Both the vertices of the octagon and the isometries are represented by complex numbers: 8(m + 1)

complex numbers for the 8(m + 1) points and 32(m + 1) complex numbers for the 8(m + 1)

isometries. Each such complex number is represented by two rational numbers and each such rational number is represented by two integers: its numerator and its denominator. The running time of a sequence of m ≥ 0 twists thus depends on the sizes of these 160(m + 1) integers; here the size of an integer is the number of digits of its decimal representation.

Twisting promptly gives rise to big numbers. As an example, take the rational admissi- Consider first the power sequence of twists 0 m (the choice of 0 is without loss of generality) for m ∈ {0, . . . , 3000}. A simple recursion gives the following for every m:

τ k,m =    (τ 0,0 ) m • τ k,0 if k ∈ {1, 2, 3} mod 8 τ k,0 • (τ 0,0 ) -m if k ∈ {5, 6, 7} mod 8 τ k,0 if k = 0 mod 4 z k,m = (τ 0,0 ) m (z k,0 ) if k ∈ {1, 2, 3} mod 8, z k,0 otherwise.
From that it is easy to see that the sizes of the integers involved in the representations of (τ k,m ) 0≤k≤7 and (z k,m ) 0≤k≤7 grow at most linearly in m. See Figure 5.

Comparison of flip strategies

As recalled in Section 2.2, a Delaunay flip algorithm can flip Delaunay-flippable edges in any order. In this section, we consider six strategies:

• naive strategy: choose the first Delaunay-flippable edge given by the cgal combinatorial map iterator DartRange::iterator.

• random strategy: choose uniformly among all the Delaunay-flippable edges.

• minimag and maximag strategies: choose the edge e whose cross-ratio R T (e) minimizes (resp. maximizes) Im [R T (e)] among the Delaunay-flippable edges.

• minratio and maxratio strategies: choose the edge e whose cross-ratio R T (e) minimizes 1. We first construct the set X containing the 11 prefixes of the sequence of twists S l (including the empty sequence) for every l ∈ {1, . . . , n s }: X contains at most 10n s + 1 sequences whose sizes vary between 0 and 10. Then for every k ∈ {1, . . . , n q } and every t ∈ X, we perform steps 2 and 4 with (z 0 , z 1 , z 2 , z 3 ) = Q k , n p interior vertices, and

(
t 1 • • • t m = t.
In each case, we run the Delaunay flip algorithm six times: one for each strategy, and count the number of flips that were needed for the algorithm to terminate. Among those six integers we denote by α k,t the minimum and by β k,t the maximum. Experiments E, F, G, and H use power sequences of twists. They are parameterized by n q , n p , and a set Ω of integers giving the lengths of the considered twists, see Table 1. For every k ∈ {1, . . . , n q }, every m ∈ Ω, and every u ∈ {0, 1, 2, 3}, we perform steps 2 and 4 with for experiments E, . . . , H) when the strategy induced the minimum/maximum number of flips among the other strategies. Figure 9 summarizes the results. Overall the minratio and the maxratio strategies seem to regularly achieve the maximum and the minimum (respectively).

(z 0 , z 1 , z 2 , z 3 ) = Q k , n p interior
Observe in particular that in experiments D and H the minratio and the maxratio strategies always induced more and fewer flips, respectively, than any other strategy. 7 Exploring the relationship between number of flips and diameter

Rationale for the experiments

Mark Bell [START_REF] Mark | Simplifying triangulations[END_REF] showed that the structure of the mapping class group has a very interesting effect on the flip graph of topological triangulations. In this topological setting, the objective is to reduce the number of intersections k of the input triangulation with a fixed curve. The main theorem of Bell's paper states that one can always find a flip or a power of a Dehn twist that reduces the number of crossings by a fixed percentage. This result can be seen as follows: either a pseudo-Anosov transformation allows the number of crossings to decrease in a single application, or there exists a power of a Dehn twist that reduces the number of crossings. This gives an algorithm to compute the optimal triangulation using O(log(k)) operations.

Our problem is different from Mark Bell's: in his study, the number of crossings is an explicit measure of the distance to the goal, while there is no way to know in advance how far the input triangulation is from being Delaunay, and we do not know the homotopy classes of final edges.

However, asymptotically, combinatorial intersection metrics are very similar to the hyperbolic metrics on surfaces of genus g ≥ 2. If a triangulation has very long edges (in terms of the number of crossings for the topological version, or in terms of the hyperbolic length in our geometric setting), then in the first stage both strategies aim at reducing edge lengths. Thus the two problems might have a similar asymptotic efficiency.

This raises two questions:

• Is there any hope to experimentally observe such similarities in the efficiency? It looks a priori unpromising as the above only holds asymptotically.

• Can Mark Bell's result be transposed to the number of flips?

We carry out two sets of experiments. The first set constructs the input triangulation by twisting the initial octagon in one direction only; as these twists correspond to reducible elements of Mod(S 2 ) (Section 2.3) we expect to observe a linear number of flips. The second set of experiments twists the octagon in a random way; asymptotically, we should obtain pseudo-Anosov elements of the mapping class group and an asymptotic logarithmic behavior.

We present five experiments named I, J, K, L, and M, all using the naive strategy (see Section 6). We use again the same notation as in Section 5. We follow steps 3 and 4 and keep track of the loosely-symmetric octagon G[z 0 , . . . , z 7 ] obtained in step 4 after the twists;

we compute (an approximation represented by a C ++ double of) the hyperbolic diameter of G[z 0 , . . . , z 7 ] (see Appendix E). As we are only interested in the influence of the diameter, we do not run step 2 (i.e., we set n p = 0) and the triangulation thus has only one vertex.

Exploring with power sequences

Experiments I and J are parameterized by the number n q of 4-tuples: n q = 1 in I and n q = 1, 000 in J. We perform step 4 with

(z 0 , z 1 , z 2 , z 3 ) = Q k , n p = 0,

Exploring with random sequences

In the following experiments the values of n q and n s are respectively n q = 1, n s = 10, 000

(experiment K), n q = 10, n s = 1, 000 (experiment L), and n q = 1, 000, n s = 100 (experiment M ). We first construct the set X containing the 11 prefixes of S l (including the empty sequence)

for every l ∈ {1, . . . , n s }. Then for every k ∈ {1, . . . , n q } and every t ∈ X, we perform step 4 with A Cross-ratios and Delaunay flips

(z 0 , z 1 , z 2 , z 3 ) = Q k , n p = 0,
We define the map φ :

C × C → C by φ(x, y) = 1 -(1 -x) • y for every (x, y) ∈ C 2 .
To prove Lemma 5, we first give a straightforward lemma. Here, the triangulation T may be infinite or finite. • R T (w) = φ(R T (w), R T (e)) for w ∈ {f, h}.

• R T (w) = φ(R T (w)), 1/R T (e )) for w ∈ {g, k}.

It is clear that the cross-ratio of any edge of T other than {e, f, g, h, k} remains unchanged after the flip.

Proof. Consider the notation defined by Figure 15. A straightforward computation gives:

[z 1 , z 2 , z 5 , z 3 ] • [z 1 , z 2 , z 3 , z 4 ] = [z 1 , z 2 , z 5 , z 4 ] [z 2 , z 3 , z 6 , z 4 ] • [z 2 , z 3 , z 4 , z 1 ] = [z 2 , z 3 , z 6 , z 1 ] [z 3 , z 4 , z 7 , z 1 ] • [z 1 , z 2 , z 3 , z 4 ] = [z 3 , z 4 , z 7 , z 2 ] [z 4 , z 1 , z 8 , z 2 ] • [z 2 , z 3 , z 4 , z 1 ] = [z 4 , z 1 , z 8 , z 3 ].
The result follows.

Let us now state the result on S.

Lemma 5. Consider a triangulation T of S and an edge e of T . Let f, g, h, k be the edges of T such that e, f, g and e, h, k (oriented counter-clockwise) bound the triangles incident to e in T .

Assume that e is Delaunay-flippable and let T be the triangulation obtained from T after the flip of e and e be the new edge replacing e. Then the following holds:

• R T (e ) = R T (e)/(R T (e) -1).

• If f = h then R T (w) = φ(R T (w), R T (e)) for every w ∈ {f, h}.

• If f = h then R T (f ) = φ(φ(R T (f ), R T (e)), R T (e)).
• If g = k then R T (w) = φ(R T (w), 1/R T (e )) for every w ∈ {g, k}.

• If g = k then R T (g) = φ(R T (g), φ(R T (g), 1/R T (e )), 1/R T (e )).
Before going to the proof, note that some edges in X = {e, f, g, h, k} may be equal. However, the edges e, f, g are pairwise-distinct and so are the edges e, h, k as they bound faces of T . Also, f = k and g = h because the interior angles of the faces of T are all less than π. Hence the only two possible equalities in X are between f and h, and between g and k.

One easily sees that the cross-ratio of any edge w / ∈ X remains unchanged after the flip.

Proof. Consider the lift T of T . Choose a fixed lift e of e and let f , g, h, k be the edges of T such that e, f , g and e, h, k bound the two faces incident to e in T , oriented counter-clockwise.

By renaming f , g to h, k and vice versa if needed we can also assume that each w ∈ X is lifted by w. We define X 1 as { f , h} if f = h or as { f } if f = h. We define X 2 similarly for g and k.

Then we set X = { e} ∪ X 1 ∪ X 2 . This way X contains exactly one lift of each element of X.

Define E as the set of all lifts of e that are incident to one of the faces of T having an edge in X.

The possible configurations are summarized in Figure 16. Consider the infinite triangulation T A 2-dimensional combinatorial map can be described as a finite set whose elements are called darts together with three permutations β 0 , β 1 , and β 2 of this set of darts. The permutations β 0 and β 1 are the inverse of each other while the permutation β 2 is an involution. We use 2dimensional combinatorial maps to describe graphs cellularly embedded on surfaces as follows.

For each face of a graph we constitute a cycle of darts such that given a dart d the next dart in the cycle is β 1 (d) (and thus the previous one is β 0 (d)). The darts of the cycle represent the edges bordering the face. We "glue" faces along their borders by pairing darts: given two darts d and d we set β 2 (d) = β 2 (d ). It is possible to identify two darts that belong to a single face.

We refer to the literature for a formal definition [START_REF] Mohar | Graphs on Surfaces[END_REF]Section 3.3]. Now we explain the role played by the anchor A in the data structure (M, F, A) described in Section 3.1. If z 1 , z 2 , z 3 , and z 4 denote 4 distinct complex numbers then the number z 4 can be deduced from z 1 , z 2 , and z 3 and from the cross-ratio [z 1 , z 2 , z 3 , z 4 ]. This fact has a consequence useful to us: given an infinite triangulation T of the hyperbolic plane if one knows the coordinates in D of the 3 vertices of some face of T together with the cross-ratio of every edge of T then one can recursively compute the coordinates of any point of T . In our setting T is the lift T of a triangulation T of a surface and for every edge e of T we have R T ( e) = R T (e) by definition. In the data structure (M, F, A) representing the triangulation T the cross-ratios of T are given by F and the anchor A provides the coordinates of the 3 vertices of some face of T . Consequently given (M, F, A) one can construct the coordinates of the vertices of T . That enables drawing some part of T for example (this is our use of the anchor).

B.2 Flipping edges (Section 3.2)

In this section, we explain how the flip of an edge e of a triangulation T is encoded on the data structure (M, F, A) defined in Section 3.1 and representing T . We may not make the distinction between the edges of the triangulation T and those of the combinatorial map M . We are given as input a dart d e of the combinatorial map M belonging to the edge e to be flipped. The edge flip is performed in three steps. Algorithms 1, 2, and 3 are performed in this order. In particular Algorithms 2 and 3 use the variables d f , d g , . . . defined in Algorithm 1. We use the notation of Section 3.1. We recall here that F denotes the function that maps each edge of the combinatorial map M to its cross-ratio and A = (δ, a 1 , a 2 , a 3 ) is the anchor. Algorithm 1 performs operations in the combinatorial map M (Figure 17), and is implemented using the cgal package [START_REF] Damiand | CGAL User and Reference Manual[END_REF].

d f ← β 1 (d e ); d g ← β 1 (d f ); d e ← β 2 (d e ); d h ← β 1 (d e ); d k ← β 1 (d h ); β 1 (d f ), β 1 (d e ), β 1 (d k ) ← d e , d k , d f ; β 1 (d g ), β 1 (d h ), β 1 (d e ) ← d h , d e , d g ; β 0 (d e ), β 0 (d k ), β 0 (d f ) ← d f , d e , d k ; β 0 (d h ), β 0 (d e ), β 0 (d g ) ← d g , d h , d e ;
Algorithm 1: Flipping the edge containing a dart d e in a combinatorial map. Recall that the lift T of T is precisely the infinite triangulation of H obtained by flipping the lifts of the edge e (there are infinitely many of them) in the lift T of T . Before the flip the anchor A represents a face t of T . We assume that t is adjacent to a lift e of the edge e to be flipped (otherwise Algorithm 2 does nothing and the anchor is correctly not modified). Let e be the edge obtained after the flip of e in T and e be the lift of e obtained after the flip of e in T . We claim that after the execution of Algorithm 2 the new anchor represents one of the two faces of T that are incident to e . First observe that such a face, name it t , shares two vertices with t and that the vertex of t that is not shared with t can be computed from the three vertices of t and the cross-ratio of e in T : this computation is done by a function φ that we now define. Let Ω ⊂ C 4 be the set of 4-tuples (x, y, z, r) ∈ C 4 such that x, y, z are pairwise-distinct and r(z -y) = (z -x).

Then φ : Ω → C is defined by φ(x, y, z, r) = (xr(z -y) + y(x -z))/(r(z -y) + x -z) on every since µ < π 6 and sin( 7π 8 ) > 1 3 . From the bounds on P 0 and P 1 , and using Assumption (2) we deduce that P 0 > P 1 , so λ is well defined. Also, we get that 1 < λ < 1 1 -3 2 µe 4R . It remains to prove that V ∈ B(f (z 3 ), ε). For the sake of clarity we denote by D the hyperbolic distance d(0 C , f (z 3 )). It is enough to show the following:

λ tanh D 2 < tanh D + ε 2 .
We first observe that x → tanh(x) -x/2 is increasing on [0, 1/2[ and maps 0 to 0. Thus, since ε/2 ∈]0, 1/2[, tanh(ε/2) ≥ ε/4. From that and by applying Assumption (2) we obtain µe 4R < tanh(ε/2)e -2R .

We conclude with the following implications: That concludes the proof.

3µe 4R <
D Details for the generation of input (Section 5)

D.1 Generating an initial rational 4-tuple (step 1)

We follow the construction of 4-tuples [1, Section 3] recalled in Section 2.4 but only compute rational approximations of the algebraic numbers involved. Then we apply Algorithm 5. The generation process described below has the following advantage: the size of the integers involved in the output 4-tuple (z 0 , z 1 , z 2 , z 3 ) are controlled by the parameter N defined below. Providing small fractions is important as the process described in this section is performed at the very beginning of the experiments.

We first construct for k ∈ {1, 2, 3} the real and imaginary parts x k and y k of a complex number z k ; they are represented as float numbers in python and constructed in [-1, 1] using the random method of the random package: x k = 2*random.random()-1. That simulates a uniform distribution. The construction fails if one of the points {z 1 , z 2 , z 3 } lies outside D, or if the condition a + b + c < 0 with a, b, and c as defined in Section 2.4 is not satisfied.

Then we construct the float numbers x 0 and y 0 representing the real and imaginary parts of z 0 as described in Section 2.4. From that we construct for each k ∈ {0, 1, 2, 3} the real and imaginary parts x k and y k of z k as rational approximations of x k and y k : we set x k = int(N * x k ) / N, where the parameter N ∈ N\{0 N } determines the quality of the approximation and int is native in python. We arbitrarily chose N = 100 in each computation. The construction fails if (z 0 , z 1 , z 2 , z 3 ) is not valid.

The rational 4-tuple (z 0 , z 1 , z 2 , z 3 ) is not necessarily admissible. However, by the construction method, it can be seen as a rational approximation of some admissible 4-tuple and it satisfies

  Poincaré disk model of H, in which H is represented as the open unit disk D of the complex plane C. Every orientation preserving isometry f : D → D can be represented by a matrix a b c d ∈ C 2×2 such that f (z) = az + b cz + d for any z ∈ D. Observe that the matrix is not unique. In addition some matrices do not represent an isometry. Given two orientation preserving isometries f and g respectively represented by matrices A and B, the product A • B represents f • g.

  of T is the (infinite) triangulation of H whose vertices and edges are the lifts of the vertices and the edges of T . A Delaunay triangulation T of S is a triangulation whose lift T is a Delaunay triangulation in H. In other words, for each face t of T and any of its lifts t, the open disk in H circumscribing t does not contain any vertex of T . Recall that circles in the Poincaré disk model correspond to circles in the complex plane.Lifting an edge e of T to some e, together with the two triangles incident to e in the lifted triangulation T , we say that e is Delaunay-flippable if the open disks of these two triangles contain the fourth vertex of the quadrilateral formed by the two triangles. In this case, the geodesic segment e that is the other diagonal of the quadrilateral is contained in it. The Delaunay flip of e in T consists in replacing e by e and projecting it back to S by ρ.Every Delaunay flip algorithm takes as input a triangulation of S and flips Delaunay-flippable edges (in any order) until there is none left. Every such algorithm terminates and outputs a Delaunay triangulation[START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF].

Figure 1 :

 1 Figure 1: (Left) A Dehn twist along the curve c modifies the blue curve as shown. (Right) A t-twist on an admissible loosely-symmetric octagon.

β 1 β 2 Figure 2 :

 22 Figure 2: A dart in a combinatorial map (bold).

  ]. Cross-ratios are suitable for a flip algorithm, due to their well-known property: assuming that the four points are oriented counterclockwise, Im [z 1 , z 2 , z 3 , z 4 ] > 0 if and only if z 4 lies inside the open disk circumscribing the triangle (z 1 , z 2 , z 3 ). Given an edge e of a triangulation T of S we consider a lift e = ( u 1 , u 3 ) of e in D and the remaining vertices u 2 and u 4 of the two faces incident to e in T , numbering vertices counterclockwise. The cross-ratio R T (e) is defined as [ u 1 , u 2 , u 3 , u 4 ]; it is independent of the choice of the lift of e, as the cross-ratio is invariant under orientation preserving isometries of D. An edge e of T is Delaunay-flippable if and only if Im [R T (e)] > 0.

4 Figure 3 :

 43 Figure 3: The triangulation T 0 and the edges e 0 , . . . , e 4 .

  For any z ∈ D and any ε > 0, B(z, ε) denotes the open ball {z ∈ D : d(z, z ) < ε} where d(•, •)

  Proof. For two reals a and b, define ]a, b[= {z ∈ C : a < Re [z] < b and Im [z] = 0}. We first choose for every k ∈ {0, 1, 2, 3} a point z k ∈ B(z k , ε) ∩ (Q + iQ), with the additional requirement that z 0 ∈]0, 1[, but without trying to satisfy the area condition A

  , the symmetric octagons of Aigon-Dupuy et al. (Section 2.4) are loosely-symmetric octagons. Identifying the opposite sides of an admissible loosely-symmetric octagon gives a surface of genus 2 [4, Theorem 1.3.5].

7 )

 7 is obtained by t-twisting (z 0 , . . . , z 7 ). For every point z in the closure of G[z 0 , . . . , z 7 ] at least z or τ t (z) lies in the closure of G[z 0 , . . . , z 7 ].

ble 4 -

 4 tuple z 0 = 10/11, z 1 = 1/2 + 1/2i, z 2 = -1/10 + 9/10i, z 3 = -3/5 + 3/5i. Then twist (z 0 , . . . , z 3 , -z 0 , . . . , -z 3 ) by m twists and get (z 0,m , . . . , z 7,m ) ∈ (Q + iQ)8 such that G[z 0,m , . . . , z 7,m ] is an admissible loosely-symmetric octagon. Denote by τ k,m the orientation preserving isometry of D mapping z k+5,m to z k,m and z k+4,m to z k+1,m , for k ∈ {0, . . . , 7}.

Figure 5 :Figure 6 :

 56 Figure 5: Size of the numerator of Re [z 1,m ] as a function of the length m of a power sequence of twists.

Figure 7 Figure 7 :

 77 Figure 7: Experiments A, B, C, and D: points (α k,t , β k,t ) for k ∈ {1, . . . , n q } and t ∈ X.

  vertices, and t 1 . . . t m = u m . Then we run the Delaunay flip algorithm for each of the six strategies. Here, the minimum and maximum number of flips are respectively denoted by α k,m,u and β k,m,u . Figure 8 shows a stronger impact of the strategy on the number of flips than experiments A, B, C, D.To compare the six strategies, we count for each experiment and each strategy the number of times (i.e., the number of pairs (α k,t , β k,t ) for experiments A, . . . , D or pairs (α k,m,u , β k,m,u )

Figure 8 :

 8 Figure 8: Experiments E, F, G, and H: points (α k,m,u , β k,m,u ), k ∈ {1, . . . , n q }, m ∈ Ω, u ∈ {0, 1, 2, 3}.

Figure 9 :

 9 Figure 9: The number of times when each strategy induced the minimum/maximum number of flips.

  flips induced by a strategy

Figure 10 :

 10 Figure 10: The number of flips induced by the naive strategy with respect to the minimum among the six strategies in experiments A, B, C, and D.

Figure 12

 12 Figure 12 illustrates a run of the program. The diameter of the initial domain is about 139 and the diameter of the final domain is smaller than 5.

  and t 1 . . . t m = u 3m for k ∈ {1, . . . , n q }, u ∈ {0, 1, 2, 3}, and m ∈ {0, . . . , 50}, and we compute the approximate hyperbolic diameter k,m,u of G[z 0 , . . . , z 7 ]. We run the Delaunay flip algorithm, counting the number α k,m,u of flips needed for the algorithm to terminate. Figure13shows the result.

Figure 13 :

 13 Figure 13: Experiments I and J: the number of flips α k,m,u with respect to the (approximate) diameter k,m,u , k ∈ {1, . . . , n q }, m ∈ {0, . . . , 50}, u ∈ {0, 1, 2, 3}.

Lemma 4 .

 4 Consider a triangulation T of H and an edge e of T . Denote by f, g, h, k the edges, oriented counter-clockwise, of the quadrilateral formed by the two triangles of T that are incident to e (hence assuming that e is incident to two bounded faces). Assume that e is Delaunay-flippable and let T be the triangulation obtained from T when replacing e by the other diagonal e of the quadrilateral. Then: • R T (e ) = R T (e)/(R T (e) -1).

Figure 15 :

 15 Figure 15: Notation for the proof of Lemma 4 (geodesic edges are represented by straight line segments).

fFigure 16 :

 16 Figure16: The possible configurations: e is the black segment, X \ { e} = X 1 ∪ X 2 is in blue, and E \ { e} is in green.of the hyperbolic plane obtained from T after flipping each element of E. We denote by e the edge of T resulting from the flip of e. Then for every w ∈ X \ { e} we have R T (w) = R T ( w) and R T (w) = R T ( w). Also, R T (e ) = R T ( e ) and R T (e) = R T ( e). The result follows by computing R T ( w) and R T ( e ) using Lemma 5.

Figure 17 :

 17 Figure 17: Illustration of Algorithm 1 In Algorithms 2 and 3 given a dart d of the combinatorial map M we denote by [d] the edge that contains d. Algorithm 2 computes an anchor for the triangulation T obtained after the flip.

2 and u 4 2 : 3 :

 2423 (x, y, z, r) ∈ Ω. The map φ is well-defined. Now we briefly explain why φ computes this third vertex of t and why Algorithm 2 always gives to φ inputs that are in Ω. Consider an infinite triangulation T of H and an edge w of T . Denote by u 1 and u 3 the vertices of w and by u the two other vertices of the two faces of T containing w: assume that u 1 , u 2 , u 3 , u 4 are in counter-clockwise order. A simple computation shows that (u 1 , u 2 , u 3 , R T (w)) ∈ Ω andu 4 = φ(u 1 , u 2 , u 3 , R T (w)).The correctness of Algorithm 2 follows by case analysis. switch δ do case d e do δ ← d h ; a 2 ← φ(a 2 , a 3 , a 1 , F ([d e ])); end case d e do δ ← d f ; a 2 ← φ(a 2 , a 3 , a 1 , F ([d e ])); end case d f or d h do a 3 ← φ(a 1 , a 2 , a 3 , F ([d e ])); end case d g or d k do a 3 ← φ(a 3 , a 1 , a 2 , F ([d e ])); end end Algorithm Updating the anchor A = (δ, a 1 , a 2 , a 3 ). The dart δ is modified if δ ∈ {d e , d f , d g , d e , d h , d k }.The update of the cross-ratios encoded in the map F is done by Algorithm 3. Algorithm 3 is a straightforward implementation of Lemma 5.F ([d f ]) ← 1 -(1 -F ([d f ])) • F ([d e ]); if β 2 (d f ) = d h then F ([d f ]) ← 1 -(1 -F ([d f ])) • F ([d e ]); else F ([d h ]) ← 1 -(1 -F (d h ])) • F ([d e ]); end F ([d e ]) ← F ([d e ])/(F ([d e ]) -1); F ([d g ]) ← 1 -(1 -F ([d g ]))/F ([d e ]); if β 2 (d g ) = d k then F ([d g ]) ← 1 -(1 -F ([d g ]))/F ([d e ]); else F ([d k ]) ← 1 -(1 -F ([d k ]))/F ([d e ]); end Algorithm Updating the cross-ratios.

Table 1 :

 1 We present eight experiments A, B, C, D, E, F, G, and H, allowing us to compare the number of flips that the six strategies induce on a variety of inputs. The notation Q k , S k and the parameters n q , n s , and n p are defined in Section 5. Parameters for experiments A to H.

	exp. A B	C	D	exp.	E	F	G	H
	n q	50 30 10	1	n q	100	30	10	10
	n s	50 30 10	10	Ω	0, 30, 60, 0, 30, 60, 0, 10,	0, 5,
						90, 120	90, 120 20, 30	10
	n p	0 10 100 1, 000	n p	0	10	100	1, 000

resp. maximizes) the quotient |Im [R T (e)]| / |R T (e)|. Let us first check that the strategy actually has an influence on the number of flips. Experiments A, B, C, and D use random sequences of twists. The values of n q , n s , and n p are shown in Table

  and t 1 . . . t m = t. We compute the approximate diameter k,t of G[z 0 , . . . , z 7 ]. We run the Delaunay flip algorithm and count the number α k,t of flips needed for the algorithm to terminate. Figure14shows α k,t as a function of 10 ln( k,t ) for k ∈ {1, . . . , n q } and t ∈ X. Here ln denotes the natural logarithm (base e). [17] Georg Osang, Mael Rouxel-Labbé, and Monique Teillaud. Generalizing CGAL periodic Delaunay triangulations. In Proceedings 28th European Symposium on Algorithms, pages 75:1-75:17, 2020. Best Paper Award (Track B: Engineering and Applications). doi:10. 4230/LIPIcs.ESA.2020.75.[18] Chee Yap et al. The CORE library project. URL: http://cs.nyu.edu/exact/core_pages/.

  tanh(ε/2)e -2R

	=⇒	3 2	µe 4R <	1 2	tanh(ε/2)(1 -tanh(R) 2 )
	=⇒	tanh(D/2) + tanh(D/2) 2 tanh(ε/2) tanh(D/2) + tanh(ε/2)	< 1 -	3 2	µe 4R
	=⇒	tanh( D 2 ) tanh( D+ε 2 )	< 1 -	3 2	µe 4R
	=⇒ λ tanh	D 2	< tanh	D + ε 2	.

Acknowledgements The authors want to thank Vincent Delecroix, Matthijs Ebbens, Hugo Parlier, Jean-Marc Schlenker, and Gert Vegter for helpful discussions over many years.

Source code available at https://members.loria.fr/Monique.Teillaud/Exp-hyperb-flips/.

Funding This work was partially supported by grant ANR-17-CE40-0033 of the French National Research Agency ANR (project SoS) https://SoS.loria.fr/.

Input:

The cross-ratios R 0 , . . . , R 4 for k = 0, . . . ,

Algorithm 4: Updating R 0 , . . . , R 4 along the sequence of 5 flips.

C.2 Approximation algorithm (Section 4.2)

This section gives additional details on the construction of the rational admissible 4-tuple shown in the proof of Theorem 2 in Section 4.2.

any k ∈ {0, 1, 2, 3} and 1 ≤ l < m ≤ 3 the following properties are satisfied:

• arg z 0 = 0

In what follows we consider some ε, µ > 0 and a rational (ε, µ)-admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ).

Algorithm 5 returns a rational admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ) such that z k ∈ B(z k , ε) for every k ∈ {0, 1, 2, 3}. We prove the correctness of the algorithm in Proposition 8 under certain assumptions on ε, µ, and on the input 4-tuple. Before describing the algorithm we state a preliminary Lemma.

Proof. This is clear since G[-z 0 , z 0 , z 2 ] is a triangle so its hyperbolic area is at most π.

Proposition 8 (Correctness of Algorithm 5). Assume ε ∈]0, 1[ and µ ∈]0, π/6[. We introduce the following parameter:

If the following assumption is satisfied:

then Algorithm 5 is well-defined and correct.

Input : Reals ε, µ > 0 and a rational (ε, µ)-admissible 4-tuple

Algorithm 5: The approximation algorithm.

Proof. We only consider the case

as the same arguments hold for the other case after application of Lemma 7. Using the notations introduced in the algorithm we bound P 0 and P 1 and then deduce that λ is well-defined. Only then we prove that V ∈ B(f (z 3 ), ε). We have

We already proved that arg

In addition,

By Definition 6, and since

Moreover for every k ∈ {0, 1, 2, 3} one has d(0

Writing the Euclidean norm in term of the hyperbolic distance

Together with the bound on arg Z 0 we obtain a bound on P 0 :

the hypothesis of Proposition 8. We can thus compute an admissible 4-tuple using Algorithm 5

(see also Section 4.2).

To simplify notation, we still denote the rational admissible 4-tuple that we obtain by

D.2 Generating points in an admissible symmetric octagon (step 2)

Consider the rational admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ) obtained after step 1. In this section, we describe our method to construct a point p ∈ Q + iQ in the closure of the admissible symmetric octagon P[z 0 , z 1 , z 2 , z 3 ], simulating a uniform distribution with respect to the hyperbolic metric.

The method uses inexact computation so it can fail especially if the Euclidean area of

This is also why we do not generate such points in the admissible loosely-symmetric octagons resulting from the twists in step 4.

We start by dividing P[z 0 , z 1 , z 2 , z 3 ] into 6 hyperbolic triangles ∆ 1 , . . . , ∆ 6 . We compute the hyperbolic area of each triangle as a C ++ (native) double number using Equality (1). Then we choose the triangle ∆ k that will contain p with probability A(∆ k )

, k ∈ {1, . . . , 6}. By a translation we can assume that 0 C is a vertex of ∆ k . We construct as double numbers the real and imaginary parts of a complex number p ∈ D, simulating a uniform choice within the closure of ∆ k . To construct p from p we cast the real and imaginary parts of p into CGAL::Gmpq numbers [START_REF] Hemmer | CGAL User and Reference Manual[END_REF]. Then we check using Lemma 9 whether p actually belongs to the closure of ∆ k ; if this is the case we return p.

and the above expression is an equality if and only if z 3 ∈ l.

Proof. The result follows from observing that τ is an orientation preserving isometry of D sending

D.3 Constructing the data structure (step 4)

After step 1, step 2, and step 3, we are given a rational admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ), points (p 1 , . . . , p np ) ∈ (Q + iQ) np lying in the closure of P[z 0 , z 1 , z 2 , z 3 ] and a sequence t 1 , . . . , t m of twists. The rational admissible 4-tuple (z 0 , z 1 , z 2 , z 3 ) defines the surface S.

Applying the procedure in Section 5.1 we construct the vertices z 0 , . . . , z 7 of the rational admissible loosely-symmetric octagon resulting from twisting P[z 0 , z 1 , z 2 , z 3 ] according to the sequence t 1 , . . . , t m . From p 1 , . . . , p np we also construct new points (p 1 , . . . , In the end we also compute the orientation preserving isometries (τ k ) 0≤k≤7 pairing the opposite sides of G[z 0 , . . . , z 7 ] (see Section 5.1).

Recursively, we construct a sequence T 0 , . . . , T np of triangulations of the octagon G[z 0 , . . . , z 7 ].

We start with the triangulation T 0 whose edges are the eight sides of G[z 0 , . . . , z 7 ], and the five geodesic segments between z 0 and z 2 , z 3 , z 4 , z 5 , z 6 . The triangulation T 0 is represented by a combinatorial map M 0 and a map P 0 associating to each vertex v of M 0 its position P 0 (v)

E Computation of the approximation of the diameter (Section 7)

Consider a rational admissible loosely-symmetric octagon O given by the 16 rational numbers representing the real and imaginary parts of its 8 vertices. The hyperbolic diameter of O is the maximum of the hyperbolic distances between any two of its vertices. For every pair z 1 , z 2 of two such distinct vertices we compute an approximation represented by a C ++ double D of the hyperbolic distance between z 1 and z 2 . The maximum (obtained using std::max) of these 8 2 values is an approximation of the hyperbolic diameter of O.

We compute every such D as follows. The isometry f : z → (z -z 1 )/(1 -z 1 z) maps Q ∩ D to a subset of Q and maps z 1 to 0. We compute the exact rational value r 2 of the square of the modulus of f (z 2 ). Then we convert r 2 to a CORE::Expr


Didier Maquin Frédéric Kratz

José Ragot

Robust Generation of Analytical Redundancy Equations and Applications to Diagnosis

faults from the model inaccuracy, we must use robust redundancy relations which are insensitive to undesired perturbations. We consider a linear system described by the following state space equation :

This paper describes a procedure for generating analytical redundancy equations for process fault detection and isolation. The proposed procedure uses a direct approach by considering the relations between the inputs and the outputs of the process ; systems with known or partially known inputs are considered and it is pointed out that the same technique yields for these different situations.

λ x(k) = A x(k) + B u(k) + B _ f(k) (1a) y(k) = C x(k) + D u(k) + D _ f(k) (1b)
where f(k) is a perturbation vector, λ states for the s or z operator depending on the representation of the model in continuous or discrete form and B _ and D _ represent the perturbation matrices of appropriate dimensions. It is assumed in the following that C is a full rank matrix. Robust redundancy equations can be generated by the elimination of the state x(k) and the perturbation vector f(k). Therefore, using (1), the equation to solve may be generally written :

Introduction

The design of fault tolerant control systems requires failures to be detected, identified and taken into account within acceptable time interval as not to affect excessively the system operation. Surveys on design methods for failure detection and isolation are given in the papers of Isermann (1984) andGertler (1991). In the past decade, interest has been focused on the use of analytic redundancy equations rather than massive redundancy. In particular, analytical redundancy equations have been designed to form and process residuals ; these residuals are closed to zero if no failure occurs and differ significantly from zero when a failure occurs.

( p T (λ) q T (λ) )             A -λ I B _ C D _ = ( 0 0) (2)
where p(λ) and q(λ) are two vectors with compatible dimensions with the state matrices. If this condition holds, multiplying equation ( 1) with (p T (λ) q T (λ)), yields the redundancy equations :

When certain inputs of the process cannot be measured, the state observation must be designed by eliminating the unknown inputs. This problem has received considerable attention in the literature. Several researchers have investigated state or input observers when some inputs are unknown ; among them, the reader is referred to Kurek (1983), El-Tohami (1983) or Miller (1984) for example.

In the study of Park (1988) a closed-loop observer that can identify simultaneously states and inputs has been developed. A systematic investigation of unknown input observers was carried out recently by Wünnenberg (1990) using the Kronecker canonical form transformation. A more simple design procedure is given by Hou (1991) using only algebraic approach. In the present work, we introduce a generation scheme of analytical redundancy equations.

(p T (λ) q T (λ)) (

        0 I y(k) -        B D u(k)) = 0 (3)
If the system (2) has no solution, due for example to an excessive number of perturbations, then one may look for an approximate solution. In order to find the optimal approximation, a certain performance index that contains a measure of the effects of the disturbances must be defined. If the perturbation matrices do not exist, the preceding formalism is still valid and allows to generate the redundancy relations.

General solution of the decoupling problem 2. Decoupling properties

To solve equation ( 2), rewrite it in the equivalent form : In practise, the process models are not perfectly known : parameter variations, unknown inputs and component faults significantly influence processes behavior. If these perturbations have not been included in the models, the model behavior differs from those of the actual process. Consequently, the residuals are as much sensitive to fault as to inherent uncertainties of the models. To separate the

(p T (λ) q T (λ)) (             A B _ C D _ -λ       I 0 0 0 ) = (0 0) (4)
Solving ( 4) corresponds to the well-known problem of the determination of the left nullspace of a binomial matrix (Gantmacher, 1977). Taking the transpose of equation ( 4 (

        A 0 C 0 B 0 D 0 -λ       I 0 0 0 )       p(λ) q(λ) =       0 0 (5)
The solution of equation ( 5) is obtained through a twostep algorithm. First, we solve (5) with D _ = 0 (i.e. D 0 = 0). Second, we transform system (5) in order to satisfy this property. Then, the first step leads to solve the following system :

Unknown inputs systems

Consider the linear system, with a p dimensional unknown input vector d(k), described by :

(λ I -A 0 ) p(λ) = C 0 q(λ) (6a) λ x(k) = A x(k) + B u(k) + G d(k) (14a) B 0 p(λ) = 0 (6b) y(k) = C x(k) + D u(k) + F d(k) (14b)
Equation (6b) implies that p(λ) belongs to the kernel of B 0 , so that the general solution of (6b) is expressed in term of an appropriate vector of constants :

where G and F are known constant matrices. Without loss of generality, we assume that :

rank ( )

G F = dim(d) (15) p(λ) = B 0 _ v(λ) (7)
The classical conditions of observability (Kurek, 1983), are not necessary because we do not try to estimate system state. With obvious definition of the matrices M, N, A _ , Q and B _ , the equations ( 14) can be expressed as :

Substituting p(λ) in (6a) and pre-multiplying first by N the left annihilator of C 0 (N C 0 = 0) and second by C 0 _ the left inverse of C 0 (C 0 _ C 0 = I) yields the two equations :

M λ x(k) + N y(k) = A _ x(k) + Q d(k) + B _ u(k) (16) N (λ I -A 0 ) B 0 _ v(λ) = 0 (8a) C 0 _ (λ I -A 0 ) B 0 _ v(λ) = q(λ) (8b)
Assuming that equation (15) holds, there exist two orthogonal matrices H and K and a non-singular matrix R such that : Equation (8a) may be solved in respect to v(λ) (Gantmacher, 1977, p. 29-30), then p(λ) is deduced from (7) and q(λ) from (8b). At the second step, consider the complete system (5) in which D 0 may always be expanded into :

Q = H         R 0 K T (17)
Introducing a new vector d _ (k) defined by :

D 0 = H         R 0 0 0 K T (9) d _ (k) = K T d(k) (18)
where H and K are two orthogonal matrices and R a regular matrix. With K T q(λ) = q _ (λ) and B _ 0 = H T B 0 , we have from the second equation of system (5) : and pre-multiplying equation ( 16) with H T , we can partitioned equation ( 16) as :

M 1 λ x(k) + N 1 y(k) = A _ 1 x(k) + R d _ (k) + B _ 1 u(k) (19a) B _ 0 p(λ) +         R 0 0 0 q _ (λ) = 0 (10) M 2 λ x(k) + N 2 y(k) = A _ 2 x(k) + B _ 2 u(k) (19b)
The first equation (19a) may be used to estimate the unknown input d _ (k) since R is regular. The second one (19b) may contain redundancies between the input u(k) and the output y(k) provide it is possible to eliminate the state x(k). For that purpose (19b) is rewritten : which may be expanded into : 11a) is then substituted in (6a) which gives :

B _ 01 p(λ) + R q _ 1 (λ) = 0 (11a) B _ 02 p(λ) = 0 (11b) (λ M 2 -A _ 2 ) x(k) = B _ 2 u(k) -N 2 y(k) (20) The solution q _ 1 (λ) of equation (
Redundancy equations are then obtained by searching the left null space of the pencil (λ M 2 -A _ 2 ) :

(λ I -A _ ) p(λ) = C 02 q _ 2 (λ) (12) p T (λ) (λ M 2 -A _ 2 ) = 0 (21) with : A _ = A 0 -C 01 R -1 B _ 01 (13)
The solution of equation ( 21) is the same as those of the equation ( 4). The redundancy equations can be described by : Thus, the remaining equations ( 11b) and ( 12) have the same structure as the system ( 6) and therefore may be solved by using the same technique. Therefore, the residual generation consists in the extraction of the p T (λ

) (B _ 2 u(k) -N 2 y(k)) = 0 (22)
The previously described method can be used to detect and localize faulty actuators. Consider a standard system where the i th input has been isolated :

u(k) -(λ-0.5) y 3 (k) = 0 2u(k) -2(1-λ) y 1 (k) -(4λ-3) y 2 (k) + (4λ-3) y 3 (k) = 0 λ x(k) = A x(k) + B i u i (k) + B _ i u _ i (k) (23a) 6. Conclusion y(k) = C x(k) + D i u i (k) + D _ i u _ i (k) (23b)
The problem of analytical redundancy equation design has been considered. We have pointed out a systematic procedure to design these equations. The technique can be used for systems with unknown inputs. where B i is the i th column of B and B _ i is the n(r-1) matrix obtained from B by deleting B i (with the same definition for D i and D _ i ). Let u i (k) be the i th entry of u(k) and u _ i (k) the (r-1) column vector obtained from u(k) by deleting u i (k). The structure of the system ( 23) is similar to the structure of the system (14). We can consider u i (k) as an unknown input and generate the redundancy relations independently of this input. If we exchange the roles of u i (k) and u _ i (k) which is therefore considered as unknown inputs, the redundancy equations depend on all but one input, therefore it make faulty actuators isolation easier. This approach should be compared with the dedicated observer approach (Frank, 1989). J. GERTLER. "Analytical redundancy in fault detection and isolation". Proc. of the IFAC / IMACS symp. on fault detection supervision and safety for technical processes, Safeprocess, p. 9-22, Baden-Baden, 1991.

Numerical example

Consider a third order system described by equation ( 14), where :

M. HOU, P.C. MULLER. "Design of robust observers for fault isolation". Proc. of the IFAC / IMACS symp. on fault detection supervision and safety for technical processes, Safeprocess, p. 295-300, Baden-Baden, 1991.

A =           1 0 0 0 0.5 0 0 0 0.5 B =           0 1 1 G =           1 1 0 R.
          1 1 0 0 1 1 0 0 1 D =           0 0 0 F =           0 0 0 J.E. KUREK. "
The state vector reconstruction for linear systems with unknown inputs". IEEE trans. on automatic control, vol. AC-28, n° 12, p. 1120AC-28, n° 12, p. -1122AC-28, n° 12, p. , 1983. Using the decomposition (17), we obtain : where a = 2 / 2 and K = 1. R.J. MILLER, R. MUKUNDAN. "On designing reduced-order observers for linear non-invariant systems subject to unknown inputs". Int. J. control, vol. 35, n°1 , p. 183-188, 1984.

H =         H 0 0 0 I with H 0 =         a -a 0 a a 0 0 0 1 J.F.
The decomposition (19) allows to define the pencil : control, vol. 48, n° 3, p. 1121-1136, 1988.

λ M 2 -A _ 2 =            
M. EL-TOHAMI, V. LOVASS-NAGY, M. MUKUNDAN. "On the design of observers for generalized state space systems using singular value decomposition". Int. J. control, vol. 38, n° 3, p. 673-683, 1983. from which we find the left orthogonal matrix :

  ) yields : redundancy relations from the system state equations. In fact, this extraction is possible by eliminating the unknown variables (state variables or perturbations). If the disturbance directions are unknown, the elimination affects only the state variables (

  References P.M. FRANK, I. WUNNENBERG. "Robust fault diagnosis using unknown input observer schemes". In Fault diagnosis in dynamic systems, R. Patton, P.M. Frank, R. Clark (eds), p. 47-97, Prentice Hall, 1989. F.R. GANTMACHER. The theory of matrices. Chelsea publishing company, 1977.

  , J.L. STEIN. "Closed-loop, state and input observer for systems with unknown inputs". Int. J.

  (1-λ) a(2λ-1.5) -a(2λ-1.5. Observer-based fault detection in dynamic systems. PhD, Duisbourg, 1990. Applying (22), the redundancy equations are then expressed :