# Computing a Dirichlet domain for a hyperbolic surface

Vincent Despré, Benedikt Kolbe, Hugo Parlier, Monique Teillaud

► To cite this version:

Vincent Despré, Benedikt Kolbe, Hugo Parlier, Monique Teillaud. Computing a Dirichlet domain for a hyperbolic surface. Journal of Computational Geometry, 2025, 15 (2), pp.40-56. 10.20382/jocg.v15i2a2 . hal-04995536

HAL Id: hal-04995536

https://hal.science/hal-04995536v1

Submitted on 18 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COMPUTING A DIRICHLET DOMAIN FOR A HYPERBOLIC SURFACE*

Vincent Despre† Benedikt Kolbe‡ Hugo Parlier§ Monique Teillaud¶

**ABSTRACT.** This paper exhibits and analyzes an algorithm that takes a given closed orientable hyperbolic surface and outputs an explicit Dirichlet domain. The input is a fundamental polygon with side pairings. While grounded in topological considerations, the algorithm makes key use of the geometry of the surface. We introduce data structures that reflect this interplay between geometry and topology and show that the algorithm runs in polynomial time, in terms of the perimeter of the initial fundamental domain and the genus of the surface.

## 1 Introduction and motivation

Hyperbolic surfaces and their moduli spaces play an ubiquitous role in mathematics, namely, through relationships with other areas including Riemannian geometry, number theory, geometric group theory and mathematical physics. Algorithms for surface groups, as combinatorial or topological objects, have a rich history dating back to Dehn. Recently, in part motivated by applications in other sciences [1, 17], there has been a push to understand hyperbolic structures on surfaces from the point of view of computational geometry.

Dealing with hyperbolic surfaces necessarily involves describing them — or even visualizing them — meaningfully. A fundamental domain (in the hyperbolic plane) with a side pairing is one way to determine a hyperbolic metric on the surface. Lengths of curves in a pants decomposition and their associated pasting parameters (so-called Fenchel-Nielsen coordinates) are another. No matter which construction or parameter set used, it is always interesting to know to which extent two different constructions output the “same” surface, where “same” can have different meanings. However, these representations, either by a fundamental domain or a set of Fenchel-Nielsen coordinates, are not unique, and determining a canonical representation is challenging for either option. In this paper, we tackle this question for fundamental domains, by computing a so-called Dirichlet domain.

<footnote>*This work was partially supported by grant ANR-17-CE40-0033 of the French National Research Agency ANR and INTER/ANR/16/11554412/SoS of the Luxembourg National Research fund FNR. https://SoS.loria.fr/</footnote>

<footnote>†Université de Lorraine, CNRS, LORIA, F-54000 Nancy, France, https://members.loria.fr/vincent.despre/</footnote>

<footnote>‡Hausdorff Center for Mathematics, University of Bonn, Germany, https://hyperbolictilings.wordpress.com/ (This work was done while this author was working at Université de Lorraine, CNRS, Inria, LORIA, France)</footnote>

<footnote>§Department of Mathematics, University of Luxembourg, https://math.uni.lu/parlier/</footnote>

<footnote>¶Université de Lorraine, CNRS, Inria, LORIA, F-54000 Nancy, France, https://members.loria.fr/monique.teillaud/</footnote>
Roughly speaking, a Dirichlet domain of a hyperbolic surface is a fundamental polygon in the hyperbolic plane, with a special point where distances to that point in the polygon correspond to distances on the surface. Another way of thinking of them is that it is a Voronoi cell associated to a lift of a single point of the surface to its universal cover $H^2$. A more precise definition is given in the next section. Note that any given surface has infinitely many Dirichlet domains up to isometry. This is namely due to the fact that (closed) hyperbolic surfaces are not homogeneous, and generically any choice for a basepoint of the Dirichlet domain will give a different domain. This is in strong contrast to, for example, flat tori. Nonetheless, when describing a surface via fundamental domains, the prize for the most relevant geometric domain undoubtedly goes to Dirichlet domains because they visualize the distance function for a given point. As far as we know, there is only one algorithm in the literature that computes a Dirichlet domain for a hyperbolic surface and a given point [19]. In the case of arithmetic surfaces, this algorithm has subsequently been improved and experiments indicate a significant speed up [18]. Unfortunately, the complexity of the algorithm is not studied and an analysis seems complicated.

The contribution of this paper is an algorithm that computes a Dirichlet domain efficiently, and its analysis. The point defining the domain is not given as input, but it is part of the output. The Dirichlet domain of a given input point can then be computed with a complexity that only depends of the genus of the surface [10]. Our main result is the following:

**Theorem 1.1.** Let $S$ be a closed orientable hyperbolic surface of genus $g$ given by a fundamental polygon of perimeter $L$ and side pairings. A Dirichlet domain for $S$ can be computed in $O((gL)^{6g-4})$ time.

A key ingredient is the use of *Delaunay triangulations* on hyperbolic surfaces, an area of research that has recently gained traction, both from an experimental and a theoretical perspective [3, 7, 4, 15, 8, 12]. Recently, it has been shown that the well-known flip algorithm that computes the Delaunay triangulation of a set of points in the Euclidean plane $E^2$ also works on a hyperbolic surface; the complexity result announced in Theorem 1.1 crucially depends on the only known upper bound on the complexity of this Delaunay flip algorithm [11]. The algorithm assumes the real RAM model; studying the algebraic numbers involved in the computations goes beyond the scope of this paper.

The paper is structured as follows: In Section 3, we give an overview of the algorithm and we present the data structure. Sections 4 and 5 explain in detail the main two steps of the algorithm. Finally, in Section 6 we conclude, using known results, the proof of Theorem 1.1 with the last step of the algorithm.

## 2 Preliminaries

We begin by recalling a collection of facts and setting notations, and we refer to [2, 5, 14] for details. The surfaces studied in this paper are assumed to be closed, orientable, and of genus $g \ge 2$. We begin with a topological surface and endow it with a hyperbolic metric to obtain a hyperbolic surface, denoted by $S$. A hyperbolic surface is locally isometric to its
universal covering space, the hyperbolic plane $\mathbb{H}^2$. Such surfaces can always be obtained by considering the quotient of $\mathbb{H}^2$ under the action of $\Gamma$, a discrete subgroup of isometries of $\mathbb{H}^2$ isomorphic to the fundamental group $\pi_1(S)$.

Let $S := \mathbb{H}^2/\Gamma$ be a closed hyperbolic surface of genus $g$ and fundamental group $\Gamma$. The projection map is denoted as $\rho: \mathbb{H}^2 \to S = \mathbb{H}^2/\Gamma$. We denote by $\tilde{x} \in \rho^{-1}(x)$ one of the lifts, to $\mathbb{H}^2$, of an object $x$ on $S$. More generally, objects in $\mathbb{H}^2$ are denoted with $\sim$.

A fundamental domain $\mathcal{F}$ for the action of $\Gamma$ is defined as a closed domain, i.e., $\mathrm{int}(\mathcal{F}) = \mathcal{F}$, such that $\Gamma\mathcal{F} = \mathbb{H}^2$ and the interiors of different copies of $\mathcal{F}$ under $\Gamma$ are disjoint. $\mathcal{F}$ contains at least one element of each orbit under the action of $\Gamma$ on $\mathbb{H}^2$, i.e., one lift of each point of $S$; its interior does not contain more than one.

For a point $\tilde{x} \in \mathbb{H}^2$, the Dirichlet domain $\mathcal{D}_{\tilde{x}}$ is defined as the Voronoi cell containing $\tilde{x}$, of the Voronoi diagram associated to the point set $\Gamma\tilde{x}$. In other words,

$$ \mathcal{D}_{\tilde{x}} = \{ \tilde{y} \in \mathbb{H}^2 \mid d_{\mathbb{H}^2}(\tilde{x}, \tilde{y}) \le d_{\mathbb{H}^2}(\tilde{x}, \Gamma\tilde{y}) \} = \{ \tilde{y} \in \mathbb{H}^2 \mid d_{\mathbb{H}^2}(\tilde{x}, \tilde{y}) \le d_{\mathbb{H}^2}(\Gamma\tilde{x}, \tilde{y}) \}, $$

where the equality is true since $\Gamma$ acts by isometries with respect to $d_{\mathbb{H}^2}$. The Dirichlet domain is a compact convex fundamental domain for $\Gamma$ with finitely many geodesic sides [2, §9.4] and is generally considered a canonical choice of fundamental domain. A property of Dirichlet domains, of interest for the conception of algorithms, is that, by the triangle inequality,

$$ \mathrm{diam}(\mathcal{D}_{\tilde{x}}) \le 2 \mathrm{diam}(S) \le 2 \mathrm{diam}(\mathcal{D}_{\tilde{x}}), $$

where $\mathrm{diam}(\cdot)$ denotes the diameter.

## 2.1 The Poincaré disk model

As they play an important part in our investigations, for completeness, we quickly recall the Poincaré disk model of the hyperbolic plane and its isometries. In this model, the hyperbolic plane $\mathbb{H}^2$ is represented as the open unit disk in the complex plane $\mathbb{C}$ and the unit circle consists of the points at infinity. The geodesics for the hyperbolic metric are either diameters of the unit disk or circle arcs meeting the unit circle orthogonally.

The representations of the fundamental group of an orientable surface we are interested in consist of orientation-preserving isometries of $\mathbb{H}^2$, represented by either of the two matrices

$$ \pm \begin{pmatrix} a & b \\ \bar{b} & \bar{a} \end{pmatrix} $$

with $a, b \in \mathbb{C}$ such that $|a|^2 - |b|^2 = 1$. The image of a point of $\mathbb{H}^2$ represented by $z \in \mathbb{C}$ is the point represented by $\frac{az+b}{bz+\bar{a}}$. Composition of isometries corresponds to multiplication of either of their representing matrices.

All isometries in the fundamental group of a surface are hyperbolic translations, which are characterized by having two distinct fixed points at infinity. The geodesic between these two points is called the geodesic axis of the translation. An isometry is a hyperbolic translation if and only if the absolute value of the trace of its matrix is larger than 2.
## 2.2 Curves, paths, and loops

A closed curve on a surface is the image of $S^1$ under a continuous map and it is said to be non-trivial (or essential) if it is not freely homotopic to a point. Similarly, a path is a continuous image of the interval $[0, 1]$, and the images of 0 and 1 are referred to as its endpoints. A loop is a path whose endpoints are equal; this endpoint is referred to as its basepoint. We refer to [14] for more details on basic topological notions concerning curves on surfaces.

For a closed curve or loop $c$, we will denote by $[c]$ its free homotopy class, and, if $c$ is a loop based in a point $p$, by $[c]_p$ its homotopy class of loops based in $p$. For a path $c$ between points $p$ and $q$, we denote by $[c]_{p,q}$ the homotopy class of the path with fixed endpoints. We will readily make use of the fact that if $c$ is a closed non-trivial curve on a hyperbolic surface, then there is a unique closed geodesic in the homotopy class $[c]$. Similarly, if $c$ is a loop based in $p$, in $[c]_p$ there is a unique geodesic loop, and if $c$ is a path between $p$ and $q$, in $[c]_{p,q}$ there is a unique geodesic path. If $c$ is a simple closed curve then the closed geodesic in $[c]$ is also simple, but this is no longer necessarily the case for loops or paths with basepoints.

The intersection number $i(c, c')$ between homotopy classes of curves $c$ and $c'$ is defined as the minimal intersection among its representatives. Note that closed geodesics on a hyperbolic surface always intersect minimally. The situation for paths is slightly different. The unique geodesic representatives of paths (with fixed endpoints) might not intersect minimally. For instance, take four points close together on a surface and such that their distance paths intersect. While the geodesics intersect it is possible to move one of them (keeping the end points fixed) by a homotopy such that they no longer intersect.

## 2.3 Fundamental polygon

Let $S$ be a (closed) hyperbolic surface of genus $g$ and fundamental group $\Gamma$. A polygon $P \subset \mathbb{H}^2$ (i.e., a circular sequence of geodesic edges) bounding a fundamental domain for $\Gamma$ (as defined in the introduction) is called a fundamental polygon. Poincaré's theorem implies that $\Gamma$ is generated by the side pairings on $P$ [2, §9.8]. The edges and vertices of $P$ project to a graph $G_P$ on $S$; the region enclosed by $P$ projects to the unique face of $G_P$.

The numbers $n_G$ of vertices and $m_G$ of edges of $G_P$ satisfy Euler's relation $n_G - m_G + 1 = 2 - 2g$, as there is only one face. It follows that if $G_P$ only has one vertex, then that vertex is incident to the $m_G = 2g$ edges, which are actually all loops. The number of vertices is maximal when they all have degree 3 (then there are no loops); in this case $3n_G = 2m_G$, so, $m_G = 6g - 3$ and $n_G = 4g - 2$. More generally, the number $2m_G$ of edges and vertices of $P$ lies between the two extreme cases: $4g \le 2m_G \le 12g - 6$. Some vertices of $P$ project to the same vertex of $G_P$, i.e., they belong to the same orbit under $\Gamma$. The graph $G_P$ has a loop for each edge of $P$ whose vertices are in the same orbit; then the projected point on $S$ is incident to that loop twice.
## 3 Algorithm overview

Let $S$ be a (closed) hyperbolic surface of genus $g$ and fundamental group $\Gamma$. We assume that, as input, we are given a fundamental polygon $\Pi \subset \mathbb{H}^2$ for $\Gamma$, together with side pairings, as in Section 2.3.

We propose the algorithm sketched below to compute a Dirichlet fundamental domain of $S$. The output of Step 1 will be denoted with primes; it will be used as input for Step 2, whose output will be denoted with double primes.

1. Find a system $\beta'_0, \dots, \beta'_{2g-1}$ of simple topological loops based at a common point $b'$ that cuts $S$ into a disk. Construct a topological polygon $\Pi'$ from lifts of these loops, together with side pairings. This step is described precisely in Section 4. The point $b'$ is given by one of its lifts $\tilde{b}'$, which can be any point in the interior of the input fundamental polygon $\Pi$.

2. Construct a point $b''$ so that the topological polygon $\Pi'$ leads to an embedded polygon $\Pi''$ when based at a lift $\tilde{b}''$ of $b''$. The construction of $b''$ is presented in Section 5.

3. Construct the Dirichlet domain of $\tilde{b}''$ (Section 6).

Obviously, the complexity of the algorithm heavily depends on the data structure used to store the objects involved in the constructions. As the algorithm actually operates in the universal covering space $\mathbb{H}^2$ of $S$, it is natural to present the data structure in $\mathbb{H}^2$. The data structure described below is actually equivalent to a combinatorial map [16, Section 3.3] on $S$, enriched with geometric information. In particular, for each vertex $x$ of $G_{\Pi}$ (the projection of $\Pi$ onto $S$, as in Section 2.3), the sequence of edges around $x$ is ordered (edges that correspond to a loop appear twice).

**Description of the input.** Notation introduced below is illustrated by Figure 1. Let a representative $\tilde{e}_i$, $i = 0, \dots, m-1$ be chosen for each couple of paired edges of $\Pi$ and denote as $\gamma_0, \dots, \gamma_{m-1} \in \Gamma$ the corresponding side pairings in $\Pi$: the other edge of the couple is $\gamma_i^{-1}\tilde{e}_i$, where $\gamma_i^{-1}$ is the inverse of $\gamma_i$. We denote the set of the $2m$ edges of $\Pi$ as $E_{\Pi}$ and the set of its $2m$ vertices as $V_{\Pi}$. We choose a representative $\tilde{v}_j$, $j = 0, \dots, n-1$ for each orbit of vertices of $\Pi$; $n$ is the number of vertices of $G_{\Pi}$.

Each element of $\Gamma$ can be represented as a word on the alphabet $\mathcal{A}_{\Gamma} = \{\mathbb{1}, \gamma_0, \dots, \gamma_{m-1}, \gamma_0^{-1}, \dots, \gamma_{m-1}^{-1}\}$, where $\mathbb{1}$ denotes the identity in $\Gamma$. Here, letters of $\mathcal{A}_{\Gamma}$ and the corresponding generators in $\Gamma$ are denoted by the same symbol; this should not cause any confusion.

The core of the data structure is a doubly linked circular list of edges of $\Pi$, which stores the combinatorial information. Additional information is necessary to store the geometry (i.e., the positions of the vertices of $\Pi$ in $\mathbb{H}^2$) and the side pairings. The data stored for each edge and vertex is constant, so the size of the data structure is $O(g)$ (we do not try to shave constants in the $O(\cdot)$).

Concretely, for each edge $\tilde{x} \in E_{\Pi}$, the data structure stores:
<figure><img src="image_2.png" /><figcaption>Figure 1: (Top) The graph $G_{\Pi}$. The arrow around vertex $v$ shows its incident edges. (Bottom) The fundamental polygon $\Pi$. Vertices $\tilde{s}, \tilde{u}, \tilde{v}$, and $\tilde{w}$ of $\Pi$ are chosen as representatives of the orbits of $s, u, v$, and $w$, respectively. The arrows show the combinatorics of the tiling $\Gamma\Pi$ at the three vertices of $\Pi$ in the orbit $\rho^{-1}v$: $\tilde{v}, \gamma_0\tilde{v} = \gamma_2\gamma_1^{-1}\tilde{v}$, and $\gamma_1^{-1}\tilde{v} = \gamma_2^{-1}\gamma_0\tilde{v}$.</figcaption></figure>
* two pointers $\text{prev}(\tilde{x})$ and $\text{next}(\tilde{x})$ that give access to the previous and next edges in $\Pi$, respectively (in counterclockwise order);

* two pointers $\text{source}(\tilde{x})$ and $\text{target}(\tilde{x})$ that give access to the source and target vertices of $\tilde{x}$ in $\Pi$, respectively (in counterclockwise order); when $\rho\tilde{x}$ is a loop in $G_{\Pi}$, $\text{source}(\tilde{x})$ and $\text{target}(\tilde{x})$ lie in the same orbit under $\Gamma$;

* a pointer to the paired edge $\text{pair}(\tilde{x})$ in $\Pi$;

* a letter $w(\tilde{x}) \in \mathcal{A}_{\Gamma}$ that encodes the relation between $\tilde{x}$ and $\text{pair}(\tilde{x})$:

$$w(\tilde{x}) = \begin{cases} 1 & \text{if } \tilde{x} = \tilde{e}_i \\ \gamma_i & \text{if } \tilde{x} = \gamma_i^{-1}\tilde{e}_i \end{cases} \quad \text{for some } i \in \{0, \dots, m-1\}.$$

$$\text{By definition, } \text{pair}(\tilde{x}) = \begin{cases} \gamma_i^{-1}\tilde{x} & \text{when } w(\tilde{x}) = 1 (\tilde{x} = \tilde{e}_i) \\ \gamma_i\tilde{x} & \text{when } w(\tilde{x}) = \gamma_i \end{cases}.$$

For each vertex $\tilde{y} \in V_{\Pi}$, the data structure stores:

* $\text{point}(\tilde{y})$, which is the representative point of its orbit: $\text{point}(\tilde{y}) = \tilde{v}_j$ for some $j \in \{0, \dots, n-1\}$;

* a word $\gamma_{\tilde{y}}$ on $\mathcal{A}_{\Gamma}$ (equivalently, $\gamma_{\tilde{y}} \in \Gamma$), which specifies the precise position $\gamma_{\tilde{y}} \text{point}(\tilde{y})$ of $\tilde{y}$ in $\mathbb{H}^2$.

The graph $G_{\Pi}$ lifts in the universal cover $\mathbb{H}^2$ to the (infinite) graph $\rho^{-1}G_{\Pi} = \Gamma\Pi$. In particular, the sequence of edges of $\Gamma\Pi$ incident to a given vertex $\tilde{v} \in \rho^{-1}v$ is a sequence of lifts of the edges incident to $v$ in $G_{\Pi}$. Each of these lifts is the image by an element of $\Gamma$ of an edge of $\Pi$. The precise positions in $\mathbb{H}^2$ of all vertices of $\Pi$ in the orbit $\rho^{-1}v$ can be computed using the information $\text{point}(\cdot)$ and $w(\cdot)$ stored in the data structure. Relations in the finitely presented group [6, Chapter 5.5] $\Gamma$ can be deduced by comparing the two sequences of edges of $\Gamma\Pi$ —clockwise and counterclockwise— around each vertex.

# 4 Constructing the initial system of simple loops

Step 1 of the algorithm described in Section 3 is reminiscent of a construction that is quite common in the topology literature: it consists in computing a spanning tree $\mathcal{T}$ of $G_{\Pi}$, then the edges of $\mathcal{T}$ are contracted, so that each vertex of $\mathcal{T}$ is merged into the root, and each edge of $G_{\Pi}$ that is not an edge of $\mathcal{T}$ is transformed into a loop based at the root.

Instead, our construction is performed in the dual of $G_{\Pi}$, which allows for a simpler presentation. Note also that topology is not enough in this work. We will actually compute the geometry of each loop by constructing lifted vertices in $\mathbb{H}^2$ and side pairings.

The main result of this section is as follows:

**Proposition 4.1.** Let $S$ be a closed orientable surface of genus $g$ and $\Pi$ a fundamental polygon of $S$ with $2m$ edges and side pairings as described in Section 3. A system of loops based at a common point on $S$, whose lifts form a topological polygon $\Pi'$ in $\mathbb{H}^2$, together with
side pairings, can be constructed in time $O(g)$. The total geodesic length of this system of loops is $O(gL)$, where $L$ denotes the perimeter of $\Pi$.

The topological polygon $\Pi'$ is given by a circular list of vertices in $\mathbb{H}^2$; its edges are homotopic with fixed endpoints to paths forming a fundamental domain. $\Pi'$ is not a fundamental polygon as in Section 2.3: the geodesic segments between its vertices (i.e., endpoints of these paths) will intersect in general, so they do not bound a fundamental domain. Section 5 will present the construction of a fundamental polygon from this topological polygon (Step 2 of the algorithm).

The algorithm that constructs $\Pi'$ proceeds in three phases:

(i). Compute a spanning tree $\mathcal{T}$ of $G_{\Pi}$.

(ii). Compute a set of loops that topologically cuts the surface into a disk.

(iii). Construct a topological polygon $\Pi'$ by lifting the loops of the previous phase. The side pairings in $\Pi'$ are also computed.

The rest of this section is devoted to proving Proposition 4.1, by detailing the construction.

*Proof.* The notation in this proof is the same as in Section 3.

Phase (i) is performed by a standard construction of a spanning tree $\mathcal{T}$ in time $O(m)$, i.e., $O(g)$. The tree has $n-1$ edges.

Let us denote by $G_{\Pi}^*$ the dual of $G_{\Pi}$. For phase (ii), we use a tree-cotree decomposition $(\mathcal{T}, \mathcal{C}, \mathcal{X})$, which is a partition of the edges of $G_{\Pi}: \mathcal{C}$ consists of the edges of $G_{\Pi} \setminus \mathcal{T}$ that are dual to the edges of a spanning tree of $G_{\Pi}^*$; $\mathcal{X}$ contains the $2g$ remaining edges. Here, $\mathcal{C}$ is empty since $G_{\Pi}$ has only one face, so $\mathcal{X} = G_{\Pi} \setminus \mathcal{T}$. Completing the $2g$ edges of $\mathcal{X}$ with the unique paths in $\mathcal{T}$ between their endpoints would lead to a basis of the homology [13]. Instead of considering $\mathcal{X}$, we actually consider the dual $\mathcal{X}^*$: since $\mathcal{C}^*$ is empty, the edges of $\mathcal{X}^*$ are already loops based at a single point and they form a basis of the homology.

Phase (iii) computes realizations in $\mathbb{H}^2$ of the $2g$ loops of $\mathcal{X}^*$. We first choose a point $\tilde{b}'$ in the interior of $\Pi$. The loops will be obtained by constructing $4g$ paths incident to $\tilde{b}'$.

Let us re-index the edges of $\Pi$ and the associated side pairings so that $\tilde{\mathcal{X}} = \{\tilde{e}_i, i = 0, \dots, 2g-1\} \cup \{\gamma_i^{-1}\tilde{e}_i, i = 0, \dots, 2g-1\}$. For each $i \in \{0, \dots, 2g-1\}$, we take a path $\tilde{\beta}'_i$ with endpoints $\tilde{b}'$ and $\gamma_i\tilde{b}';$ the path $\gamma_i^{-1}\tilde{\beta}'_i$ has endpoints $\tilde{b}'$ and $\gamma_i^{-1}\tilde{b}'$. The two paths $\tilde{\beta}'_i$ and $\gamma_i^{-1}\tilde{\beta}'_i$ project onto $S$ on the same loop $\beta'_i$ based at $b'$, and which lies in the homotopy class of the dual edge $e_i^*$ of $e_i \in \mathcal{X}$ (Figure 2).

The order in which the topological loops $\beta'_0, \dots, \beta'_{2g-1}$ appear around $b'$ on $S$ is simply given by the order of the elements of $\tilde{\mathcal{X}}$ around $\Pi$, which is obtained in $O(g)$ operations by traversing the boundary of $\Pi$ (as mentioned in Section 2.3, each loop appears twice in the sequence). However, in general the geodesic segments between the endpoints of each path may intersect and not appear in the right order around $\tilde{b}'$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Edges of Π, which are geodesics in $H^2$, are represented as straight lines. Lifts of $e_0$ and $e_j$, which belong to $\mathcal{X}$, are shown in green, while edges of Π that belong to $\tilde{T}$ are dashed. Lifts incident to $\tilde{b}'$ of the loops $\beta_0, \dots, \beta_{2g-1}$ corresponding to edges of $\tilde{\mathcal{X}}$ appear as blue curves.</figcaption></figure>

To construct $\Pi'$, we will walk in the faces of the graph $\Gamma\Pi$ and construct $4g$ lifts of $\beta'_0, \dots, \beta'_{2g-1}$ and their vertices in the orbit $\Gamma b'$, which will form the ordered sequence of edges of the topological polygon $\Pi'$. Along the construction, we will create a table $t$ that stores for each vertex $\tilde{x}$ of $\Pi'$ the element $\gamma \in \Gamma$ (as a word) such that $\tilde{x} = \gamma \tilde{b}$. We describe the walk in the next paragraph.

Without loss of generality, we start with $\tilde{\beta}'_0$ as a first edge of $\Pi'$, incident to $\tilde{b}'$. Its other vertex is $\gamma_0 \tilde{b}'$; we store 1 and $\gamma_0$ in $t$ for the first two vertices of $\Pi'$. The next edge of $\Pi'$ is incident to the endpoint $\gamma_0 \tilde{b}'$ of $\tilde{\beta}'_0$. It is given by the dual of the first edge of $\Pi$—when following $\Pi$ in counterclockwise order—that is an element of $\tilde{\mathcal{X}}$, which is found in the data structure presented in Section 3 as next(next(...next(pair(\tilde{e}_0))))...) (see Figure 3). Here the iteration on the edges of $\Pi$ stops as soon as it encounters an element of $\tilde{\mathcal{X}}$; it is an edge $\tilde{e}_j$ or $\gamma_j^{-1} \tilde{e}_j$, for some $j \in 0, \dots, 2g-1$. The edge of $\Pi'$ following $\tilde{\beta}'_0$ after $\gamma_0 \tilde{b}'$ is thus $\gamma_0 \tilde{\beta}'_j$ or $\gamma_0 \gamma_j^{-1} \tilde{\beta}'_j$ and its next vertex is $\gamma_0 \gamma_j \tilde{b}'$ or $\gamma_0 \gamma_j^{-1} \tilde{b}'$. We store either $\gamma_0 \gamma_j$ or $\gamma_0 \gamma_j^{-1}$ in $t$. We repeat the same operations for each edge of $\Pi'$. The total number of accesses to next(·) during the process is exactly the number of edges of $\Pi$, and there are a constant number of other operations for each edge, so, the boundary of $\Pi'$ is constructed in $O(g)$ time.

The side pairing between any two paired sides is easily obtained by comparing the words in $t$ that yield their respective vertices when applied to $\tilde{b}'$.

The homotopy class of each loop $\tilde{\beta}'_i$ contains the path formed by the projection onto $S$ of the path $\tilde{e}'_i$ formed by the shortest path in $\Pi$ between $\tilde{b}'$ and a point on $\tilde{e}_i$ (e.g., the midpoint), followed by the shortest path in $\gamma_i\Pi$ between this (mid)point and the other vertex $\gamma_i\tilde{b}'$ of $\tilde{\beta}'_i$ (see Figure 4). These two shortest paths are contained respectively in $\Pi$ and $\gamma_i\Pi$, so, each of them is not longer than the diameter of $\Pi$, which is bounded by half the perimeter $L$ of $\Pi$. Thus the geodesic length of $\tilde{\beta}'_i$, for $i = 0, \dots, 2g-1$, is not larger than $L$. $\square$
<figure><img src="image_4.png" /><figcaption>Figure 3: Construction of the edges of $\Pi'$, colored in orange.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Bounding the geodesic length of $\beta'_j$.</figcaption></figure>
*Remark 4.2.* The topological polygon $\Pi'$ is obtained as a circular sequence of images of some lifts $\tilde{\beta}'_i, i = 0, \dots, 2g-1$ by elements of $\Gamma$. The paths $c'_i$ (see the last paragraph of the proof) only intersect edges of $\mathcal{X}$, they do not intersect any edge of $\mathcal{T}$. The shape formed by the sequence of the images of all shortest paths $\tilde{c}'_i$ by the same elements of $\Gamma$ contains a lift of the tree $\mathcal{T}$ that is a tree in $\mathbb{H}^2$.

It is crucial to notice that the point $\tilde{b}'$ from which the topological polygon $\Pi'$ is constructed can be any point in the interior of $\Pi$. The following section shows how to choose a point that yields a polygon embedded in $\mathbb{H}^2$.

# 5 Finding an embedded system of loops

We want to find a collection of geodesic loops on a hyperbolic surface $S$, all based in a single point and disjoint otherwise, such that the complementary region of the loops lifts as a convex hyperbolic polygon. What we show is that in fact we can retain the choice of topological loops $\beta'_0, \dots, \beta'_{2g-1}$ made in Section 4 by choosing the basepoint appropriately to ensure that their geodesic realizations satisfy the desired properties.

Consider the set of topological loops $\beta'_0, \dots, \beta'_{2g-1}$ all based at point $b'$ constructed in the previous section. We choose a pair that intersects minimally exactly once, which, up to reordering, we suppose are $\beta'_0$ and $\beta'_1$. For future reference we set $L_0 := \max\{\ell(\beta'_0), \ell(\beta'_1)\}$, where $\ell$ denotes the length.

*Remark 5.1.* While all of the loops are based at the same point, they each correspond to free homotopy classes of closed curves that may or may not intersect. We can actually fix any loop to be $\beta'_0$ and find a loop $\beta'_1$ (essentially) intersecting it exactly once (meaning that the corresponding closed geodesic curves intersect exactly once). Indeed, the set $\beta'_0, \dots, \beta'_{2g-1}$ contains curves that pairwise intersect at most once, and are all non-separating and thus homologically non-trivial. As it generates homotopy, it also generates homology and in particular every curve must be intersected by at least one other curve. As they can intersect at most once, they intersect exactly once.

We begin by taking the unique geodesic loops, based in $b'$, in the homotopy classes of $\beta'_0$ and $\beta'_1$, and we replace the curves with these geodesic representatives (we keep the same notation for convenience). Now we further consider the unique simple closed geodesic representatives in the free homotopy class of $\beta'_0$ and $\beta'_1$, which we denote $\beta_0''$ and $\beta_1''$, respectively. By hypothesis, they intersect in a single point $b''$, which will be our new basepoint.

We now define a path between $b''$ and $b'$ as follows. We consider a single lift $\tilde{\beta}'_0$ of $\beta'_0$. Its endpoints both correspond to distinct lifts of $b'$ which are related by the translation $\gamma_0$ in $\Gamma$. The copies of $\tilde{\beta}'_0$ by iterates of $\gamma_0$ form a piecewise geodesic line $\tilde{\beta}'_0$ with the same endpoints at infinity as the geodesic axis of $\gamma_0$. This line $\tilde{\beta}'_0$ separates $\mathbb{H}^2$ into two half-spaces, only one of which is convex —unless $\tilde{\beta}'_0$ is smooth at $b'$, in which case the two half-spaces are convex. We now choose an endpoint of $\tilde{\beta}'_0$ and consider the lift of $\beta'_1$ that lies in the convex half-space and shares the same endpoint; in the case when the two half-spaces are convex, we can take any of them. This lift we denote by $\tilde{\beta}'_1$ and, as before, we consider the corresponding translation $\gamma_1$ or $\gamma_1^{-1}$ (which we assume to be $\gamma_1$, without loss of generality)
in $\Gamma$, its geodesic axis, and its corresponding piecewise geodesic line $\hat{\beta}'_1$. Now, we obtain $\tilde{b}''$ as the intersection of the axes of $\gamma_0$ and $\gamma_1$. We consider the unique geodesic path $\tilde{c}$ between $\tilde{b}'$ and $\tilde{b}''$ and its projection $c$ on $S$.

<figure><img src="image_6.png" /><figcaption>Figure 5: $c$ is homotopic on $S$ to the projection of the concatenation of the bold arc of $\tilde{\beta}'_0$ from $\tilde{b}'$ and the bold segment of the axis of $\gamma_1$.</figcaption></figure>

We observe that the axis of $\gamma_0$ must lie in an $R$ neighborhood of $\tilde{\beta}'_0$ where $R < \ell(\beta'_1)$. In particular, the axis of $\gamma_1$ intersects $\tilde{\beta}'_0$ (see Figure 5). Similarly, the axis of $\gamma_0$ intersects $\tilde{\beta}'_1$. By following an arc of $\tilde{\beta}'_0$ from $\tilde{b}'$ and then a segment on the axis of $\gamma_1$, we obtain a path between $\tilde{b}'$ to $\tilde{b}''$.

Observe that for $i = 0, 1, \beta''_i$, based in $b''$, is freely homotopic to $c^{-1} \cdot \beta'_i \cdot c$ and that there is a homeomorphism of $S$, isotopic to the identity, that takes $b'$ to $b''$ and that sends (the homotopy class of) $\beta'_i$ to $\beta''_i$. This homeomorphism is often referred to as the point pushing map (see for instance [14, Section 4.2] for details).

We can apply this same homeomorphism to the remaining curves. For $i = 0, \dots, 2g-1$ we set the homotopy class of loop $\beta''_i$ to be:

$$[\beta''_i]_{b''} = [c^{-1} \cdot \beta'_i \cdot c]_{b''}. \qquad (5.1)$$

As we have just moved the basepoint by a homeomorphism, the homotopy classes $[\beta''_i]_{b''}$ all have simple representatives and can be realized disjointly outside of $b''$. The following lemma implies that their unique geodesic representatives enjoy this same property. It is well known to specialists, but we include a proof sketch for completeness.

**Lemma 5.2.** Let $\Sigma$ be a hyperbolic surface with piecewise-geodesic boundary such that the interior angles on the singular points $s_0, \dots, s_{k-1}$ of the boundary are cone points of angle $\le \pi$. If $[\alpha]_{p_i, q_i}, [\alpha']_{p_j, q_j}$ are simple homotopy classes of paths (with endpoints $p_i, p_j, q_i, q_j$
in the set $s_0, \dots, s_{k-1}$), and disjoint except for possibly in their endpoints, then the unique geodesic representatives are also simple and disjoint.

**Sketch of proof.** We consider $\tilde{\Sigma}$, the universal cover of $\Sigma$, which we view as a (geodesically convex) subset of $\mathbb{H}^2$. We lift $\partial\Sigma$ to $\tilde{\Sigma}$ and representatives of $[\alpha]_{p_i,q_i}$ and $[\alpha']_{p_j,q_j}$, which are simple and disjoint, to the universal cover. Observe that being simple and disjoint is equivalent to all individual lifts in $\mathbb{H}^2$ being simple and pairwise disjoint. Now take two individual lifts of either $\alpha$ or $\alpha'$, and their unique geodesic representatives. We will see that they are also disjoint. Note that in general, given two simple disjoint paths in the hyperbolic (or Euclidean) plane, the unique geodesics between their endpoints might intersect (as already mentioned in Section 2.2). However:

**Observation:** Let $C \subset \mathbb{H}^2$ be a convex with non-empty boundary, and $p_0, q_0, p_1, q_1 \in \partial C$. Let $\alpha_1 : [0,1] \to C$ and $\alpha_2 : [0,1] \to C$ be simple paths, disjoint in their interior, with $\alpha_0(0) = p_0, \alpha_0(1) = q_0$ and $\alpha_1(0) = p_1, \alpha_1(1) = q_1$. Then the unique geodesic between $p_0$ and $q_0$ and the unique geodesic between $p_1$ and $q_1$ are disjoint in their interior as well.

A key point is that, thanks to the angle condition on the cone points, $\tilde{\Sigma}$ is a convex region of $\mathbb{H}^2$. (This is just a slightly more sophisticated observation than the elementary fact that a polygon with all interior angles less than $\pi$ is convex.) The observation now implies that the lifts of geodesics corresponding to $\alpha$ and $\alpha'$ are disjoint in their interior if and only if there are representatives of $[\alpha]_{p_i,q_i}$ and $[\alpha']_{p_j,q_j}$ that are, too, which, by hypothesis, is the case. $\square$

We can now apply Lemma 5.2 to the geodesic representatives of $[\beta_i'']_{b''}$. For simplicity we denote by $\beta_i''$ the unique geodesic loop in the corresponding homotopy class.

**Theorem 5.3.** Let $\beta_0'', \dots, \beta_{2g-1}''$ be a set of topological loops based in $b''$ that cuts a surface $S$ into a disk. Assume that $\beta_0''$ and $\beta_1''$ are closed geodesics. Then, the geodesic loops homotopic to $\beta_0'', \dots, \beta_{2g-1}''$ are simple and pairwise disjoint in their interiors. Furthermore, by cutting $S$ along those geodesics and lifting to $\mathbb{H}^2$, one obtains a convex hyperbolic polygon with 4g edges.

*Proof.* As $\beta_0''$ and $\beta_1''$ are closed geodesics, they form 4 angles in $b''$, and the opposite ones are equal. These angles thus satisfy $2\theta + 2\theta' = 2\pi$ so in particular both $\theta$ and $\theta'$ are strictly less than $\pi$. Thus by cutting along $\beta_0''$ and $\beta_1''$, we obtain a genus $g-1$ surface with a boundary consisting of 4 geodesic segments, and with 4 cone point singularities of angles $< \pi$ (see Figure 6).

We now proceed inductively for $i \ge 2$ and consider the unique geodesic path $\beta_i''$, which by virtue of Lemma 5.2, has disjoint interior from the previous geodesic segments. Furthermore, as each segment further splits the angles, the angles are all less than $\pi$.

The end result is a polygon with all interior angles less than $\pi$ which, by elementary hyperbolic geometry, is convex. $\square$

**Proposition 5.4.** Let $S$ be hyperbolic of genus $g$ and $\Pi'$ a topological fundamental polygon of $S$ with 4g edges and side pairings as described at the end of Section 4. A convex fundamental polygon $\Pi''$ with its side pairing and whose vertices project to a single vertex on $S$, can be
<figure><img src="image_7.png" /><figcaption>Figure 6: A visualization of the cutting along $\beta_0''$ and $\beta_1''$.</figcaption></figure>

constructed in $O(g)$ time. The perimeter of $\Pi''$ is $O(gL)$, where $L$ denotes the perimeter of $\Pi$.

*Proof.* We need to compute the output convex polygon $\Pi''$ i.e., $4g$ lifts of $b''$ and $2g$ side pairings $\gamma_0''$, \dots, $\gamma_{2g-1}''$. As homotopy classes of $\beta_i''$ and $\beta_i'$ are conjugates for $i = 0, \dots, 2g-1$ (Equality 5.1), the side pairing $\gamma_i''$ is equal to $\gamma_i'$ for each $i$.

The key point here is the computation of a lift of $b''$. The first step consists in finding the loops $\beta_0'$ and $\beta_1'$ satisfying $i(\beta_0', \beta_1') = 1$. As shown in Remark 5.1, we can choose any loop for $\beta_0'$. We also fix $\tilde{b}_0'$ to be an endpoint of one of the two paired sides of $\Pi'$ that are lifts of $\beta_0'$. We know the ordered sequence of loops around $b'$ (see the proof of Proposition 4.1); recall that each loop $\beta_0', \dots, \beta_{2g-1}'$ appears twice in the sequence (Section 2.3). We take as $\beta_1'$ one of the loops that alternate with $\beta_0'$ in the sequence, and choose for $\tilde{\beta}_1'$ one of its two lifts that are incident to $\tilde{b}_0'$.

The second step consists in finding the free geodesics in the homotopy classes of $\beta_0'$ and $\beta_1'$, respectively. Using $t$, we find the word $\gamma_0$ on $\{\gamma_0', \dots, \gamma_{2g-1}'\}$ representing the translation that sends $\tilde{b}_0'$ to the other endpoint of $\tilde{\beta}_0'$ (see Figure 7). The sequences $\gamma_0^n(\tilde{b}_0')$ and $\gamma_0^{-n}(\tilde{b}_0')$ converge in $\mathbb{C}$ to two points on the unit circle: these points are the two (infinite in $\mathbb{H}^2$) fixed points of the translation $\gamma_0$, i.e., the two solutions of equation $\gamma_0(z) = z$ in $\mathbb{C}$. The axis of $\gamma_0$, i.e., the geodesic between these two points, projects onto $S$ to the free geodesic in $[\beta_0']$.

We repeat the same process with $\tilde{\beta}_1'$ and find the geodesic in $\mathbb{H}^2$ that projects to the free geodesic in $[\beta_1']$. The point $\tilde{b}_0''$ comes as the intersection point of the two geodesics in $\mathbb{H}^2$. We now define $\tilde{\beta}_0''$ as the geodesic segment between $\tilde{b}_0''$ and $\gamma_0(\tilde{b}_0'')$, and $\tilde{\beta}_1''$ in the same way. This step is performed in constant time.

We can now compute the $4g$ lifts of $b''$ that are the vertices of $\Pi''$ by applying the elements of $t$ to $\tilde{b}_0''$. This last step has complexity $O(g)$.

Recall that the upper bound on the geodesic perimeter of $\Pi'$ given by Proposition 4.1 holds for the geodesic perimeter of $\Pi'$ constructed from any basepoint $b'$ in the interior of $\Pi$. So, it also holds for the new basepoint $b''$ and gives a bound on the perimeter of $\Pi''$. $\square$
<figure><img src="image_8.png" /><figcaption>Figure 7: The computation of $b''$ from $b'$.</figcaption></figure>

# 6 Finding a Dirichlet domain from an embedded system of loops

We first summarize what we have obtained so far. We started with a polygon $\Pi$ of perimeter $L$ and we obtained a convex polygon $\Pi''$ of total length $O(gL)$. Additionally, all vertices of $\Pi''$ project on a single vertex $b''$ on $S$. This construction has complexity $O(g)$ by Propositions 4.1 and 5.4.

For the sake of clarity, the two steps 1 and 2 of the algorithm (Section 3) were distinct and were presented in different sections (respectively, Sections 4 and 5). However, we can remark that they could have been merged by placing the point $b'$, which is left undefined in Section 4, directly at $b''$ as defined in Section 5. This is probably what would be done in practice.

Now it is easy to compute a Dirichlet domain. Indeed, we can triangulate $\Pi''$ easily since it is convex and, thus, we obtain a geometric triangulation $T$, on to which the Delaunay flip algorithm can be applied [11]. The complexity of this algorithm depends on the diameter of $T$, for which the perimeter of $\Pi''$ is an upper bound.

*Remark 6.1.* *Remark 4.2* can be expressed in a much simpler way at this point: $\Pi''$ itself contains a lift $\tilde{T}$ of $T$. As $\Pi''$ is convex, the geodesic segments dual to the edges of $\tilde{T}$ are contained in $\Pi''$ and induce a partition of $\Pi''$. We can actually use these segments to triangulate $\Pi''$. In this way, in the case when the input polygon $\Pi$ is already a Dirichlet domain, we already have a Delaunay tesselation, and no flip is performed by the algorithm.

The output of the flip algorithm is a Delaunay triangulation $DT$ of $S$ with the single vertex $b''$ computed in Section 5. To obtain a Dirichlet domain from $DT$, we just have to compute the triangles of $\overline{DT}$ incident to a lift $\tilde{b}''$ of $b''$ and their dual: we compute the circumcenter of each triangle to get the vertices of the Dirichlet domain and we put a
geodesic between vertices that correspond to adjacent triangles around $\tilde{b}''$. This step is also
clearly done in $O(g)$ operations. Putting all together we obtain the following theorem:

**Theorem 6.2.** Let $S$ be a closed orientable hyperbolic surface of genus $g$ given by a fundamental polygon of perimeter $L$ and side pairings. A Dirichlet domain of $S$ can be computed in time $O(f(gL) + g)$ where $f(\Delta)$ is the complexity of the flip algorithm for a triangulation of diameter $\Delta$ with a single vertex.

Using the best known bound $O(\Delta^{6g-4})$ for the flip algorithm so far [11], we obtain
Theorem 1.1 stated in the introduction as a corollary. Note that the constant in the $O(\cdot)$
depends on the metric on $S$. However, there are experimental and theoretical insights
suggesting that the actual complexity of the flip algorithm may be much better [9].

**Acknowledgements.** The authors wish to thank the SoCG anonymous reviewer who suggested this simpler version of Section 4.

References

[1] N.L. Balazs and A. Voros. Chaos on the pseudosphere. *Physics Reports*, 143(3):109-240, 1986. doi:10.1016/0370-1573(86)90159-6.

[2] Alan F. Beardon. *The Geometry of Discrete Groups*. Springer-Verlag, 1983.

[3] Mikhail Bogdanov, Olivier Devillers, and Monique Teillaud. Hyperbolic Delaunay complexes and Voronoi diagrams made practical. *Journal of Computational Geometry*, 5:56-85, 2014. doi:10.20382/jocg.v5i1a4.

[4] Mikhail Bogdanov, Monique Teillaud, and Gert Vegter. Delaunay triangulations on orientable surfaces of low genus. In Sándor Fekete and Anna Lubiw, editors, *32nd International Symposium on Computational Geometry (SoCG 2016)*, volume 51 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 20:1-20:17, Dagstuhl, Germany, 2016. Schloss Dagstuhl - Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.SoCG.2016.20.

[5] Peter Buser. *Geometry and spectra of compact Riemann surfaces*. Modern Birkhäuser classics. Birkhäuser, Boston, Mass., 2nd edition, 2010.

[6] H. S. M. Coxeter and W. O. J. Moser. *Generators and Relations for Discrete Groups*. Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1957.

[7] Jason DeBlois. The centered dual and the maximal injectivity radius of hyperbolic surfaces. *Geometry and Topology*, 19(2):953–1014, 2015. doi:10.2140/gt.2015.19.953.

[8] Jason DeBlois. The Delaunay tessellation in hyperbolic space. *Mathematical Proceedings of the Cambridge Philosophical Society*, 164(1):15–46, 2018. doi:10.1017/S0305004116000827.
[9] Vincent Despré, Loïc Dubois, Benedikt Kolbe, and Monique Teillaud. Experimental analysis of Delaunay flip algorithms on genus two hyperbolic surfaces. Preprint, INRIA, May 2021. URL: https://hal.inria.fr/hal-03462834/.

[10] Vincent Despré, Benedikt Kolbe, and Monique Teillaud. Representing infinite hyperbolic periodic Delaunay triangulations using finitely many Dirichlet domains. *Discrete & Computational Geometry*, 72:1–28, 2024. doi:10.1007/s00454-024-00653-x.

[11] Vincent Despré, Jean-Marc Schlenker, and Monique Teillaud. Flipping geometric triangulations on hyperbolic surfaces. In Sergio Cabello and Danny Z. Chen, editors, *36th International Symposium on Computational Geometry (SoCG 2020)*, volume 164 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 35:1–35:16, Dagstuhl, Germany, 2020. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.SoCG.2020.35.

[12] Matthijs Ebbens, Hugo Parlier, and Gert Vegter. Minimal Delaunay triangulations of hyperbolic surfaces. In Kevin Buchin and Éric Colin de Verdière, editors, *37th International Symposium on Computational Geometry (SoCG 2021)*, volume 189 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 31:1–31:16, Dagstuhl, Germany, 2021. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.SoCG.2021.31.

[13] David Eppstein. Dynamic generators of topologically embedded graphs. In *Proceedings of the fourteenth annual ACM-SIAM symposium on Discrete algorithms*, pages 599–608, 2003. URL: https://dl.acm.org/doi/10.5555/644108.644208.

[14] Benson Farb and Dan Margalit. *A Primer on Mapping Class Groups (PMS-49)*. Princeton University Press, 2012. URL: http://www.jstor.org/stable/j.ctt7rkjw.

[15] Iordan Iordanov and Monique Teillaud. Implementing Delaunay triangulations of the Bolza surface. In Boris Aronov and Matthew J. Katz, editors, *33rd International Symposium on Computational Geometry (SoCG 2017)*, volume 77 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 44:1–44:15, Dagstuhl, Germany, 2017. Schloss Dagstuhl–Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.SoCG.2017.44.

[16] Bojan Mohar and Carsten Thomassen. *Graphs on Surfaces*. Johns Hopkins University Press, Baltimore, 2001.

[17] Nikolai C Passler, Xiang Ni, Guangwei Hu, Joseph R Matson, Giulia Carini, Martin Wolf, Mathias Schubert, Andrea Alù, Joshua D Caldwell, Thomas G Folland, et al. Hyperbolic shear polaritons in low-symmetry crystals. *Nature*, 602(7898):595–600, 2022. doi:10.1038/s41586-021-04328-y.

[18] James Rickards. Improved computation of fundamental domains for arithmetic Fuchsian groups. *Mathematics of Computation*, 91:2929–2954, 2022. doi:10.1090/mcom/3777.

[19] John Voight. Computing fundamental domains for Fuchsian groups. *Journal de Théorie des Nombres de Bordeaux*, 21(2):467–489, 2009. doi:10.5802/jtnb.683.