Jean-Daniel Boissonnat

Siargey Kachanovich

{1

Tracing Isomanifolds in R d in Time Polynomial in d using Coxeter-Freudenthal-Kuhn Triangulations

Keywords: Subject Classification Theory of computation → Computational geometry Coxeter triangulation, Kuhn triangulation, permutahedron, PL-approximations, isomanifolds/solution manifolds/isosurfacing { 1, 2 }, { 3 } {2}, {1, 3} {2, 3}, {1} { 3 }, { 1, 2 } {1, 3}, {2} {1}, {2, 3} { 1, 2 }, { 3 } {2}, {1, 3} {2, 3}, {1} { 3 }, { 1, 2 } {1, 3}, {2} {1}, {2, 3}

Isomanifolds are the generalization of isosurfaces to arbitrary dimension and codimension, i.e. submanifolds of R d defined as the zero set of some multivariate multivalued smooth function

, where n is the intrinsic dimension of the manifold. A natural way to approximate a smooth isomanifold M is to consider its Piecewise-Linear (PL) approximation M based on a triangulation T of the ambient space R d . In this paper, we describe a simple algorithm to trace isomanifolds from a given starting point. The algorithm works for arbitrary dimensions n and d, and any precision D. Our main result is that, when f (or M) has bounded complexity, the complexity of the algorithm is polynomial in d and δ = 1/D (and unavoidably exponential in n). Since it is known that for δ = Ω(d 2.5 ), M is O(D 2 )-close and isotopic to M, our algorithm produces a faithful PL-approximation of isomanifolds of bounded complexity in time polynomial in d. Combining this algorithm with dimensionality reduction techniques, the dependency on d in the size of M can be completely removed with high probability. We also show that the algorithm can handle isomanifolds with boundary and, more generally, isostratifolds. The algorithm for isomanifolds with boundary has been implemented and experimental results are reported, showing that it is practical and can handle cases that are far ahead of the state-of-the-art.

Introduction

Given a surface represented in R 3 as the zero set of a function f : R 3 → R, the goal of isosurfacing is to find a piecewise linear (PL) approximation of the surface. This question naturally extends to isomanifolds of higher dimensions and codimensions defined as the zero set of multivariate multivalued smooth functions f : R d → R d-n . Isosurfaces play a crucial role in medical imaging, computer graphics and geometry processing [START_REF] Newman | A survey of the marching cubes algorithm[END_REF]. Higher dimensional isomanifolds are also of fundamental importance in many fields like statistics [START_REF] Chen | Solution manifold and its statistical applications[END_REF], dynamical systems [START_REF] Michael | The computation of fixed points and applications[END_REF], econometrics, or mechanics [START_REF] Newman | A survey of the marching cubes algorithm[END_REF].

State-of-the-art. The most widely used algorithm to trace isomanifolds is the Marching Cube (MC) algorithm and its numerous variants [START_REF] Gomes | Implicit Curves and Surfaces: Mathematics, Data Structures and Algorithms[END_REF][START_REF] Wenger | Isosurfaces: geometry, topology, and algorithms[END_REF]. The MC algorithm uses a cubical grid to tesselate the ambient space. In many applications in 3-dimensions, the ambient space is decomposed into unstructured tetrahedral meshes, which led to the development of a variant of the MC algorithm named the Marching Tetrahedra algorithm. In higher dimensions, any tessellation of the ambient space has a complexity that depends exponentially on the ambient dimension. Hence a key to extending marching algorithms to higher dimensions is to circumvent the curse of dimensionality by using an implicit representation of the ambient tessellation. This is impossible for general triangulations but easy to do if one uses a grid.

However, using a grid has other drawbacks and is not sufficient to break the exponential barrier. The reason for this is that the number of configurations inside a cubical cell grows exponentially with the dimension [START_REF] Wenger | Isosurfaces: geometry, topology, and algorithms[END_REF]. Hence the most promising approach seems to be to subdivide the ambient space R d using a highly regular triangulation such as the Freudenthal-Kuhn triangulation. Some early work along this direction has been published in Applied

Mathematics [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF][START_REF] Eaves | A course in triangulations for solving equations with deformations[END_REF][START_REF] Michael | The computation of fixed points and applications[END_REF], and a slightly more recent paper by Dobkin et al. [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF] attracted the interest of the Computer Graphics community to the related Coxeter triangulations. Dobkin et al. however only considered the case of curves (n = 1). The most advanced work we are aware of is due to Min [START_REF] Min | Simplicial isosurfacing in arbitrary dimension and codimension[END_REF]. Min's method uses the Freudenthal-Kuhn triangulation over a dyadic grid of R d and applies to isomanifolds of any dimension and codimension. The time complexity of Min's method is, with our notations, O(δ n log δ), where δ = 1/D and D is the maximal diameter of the simplices. The ambient dimension d is a constant hidden in the big O. The fact that the exponent of δ is the intrinsic dimension n, and not the ambient dimension d is a clear improvement over earlier methods. However, although not explicitly analysed by Min, the complexity in d remains exponential, and the method seems to be limited to small ambient dimensions. Experimental results are only reported in 3, and 4D.

Background

Permutahedral representation of CFK-triangulations

In this section, we give the most important definitions and basic properties of Coxeter and

Freudenthal-Kuhn triangulations. An extensive discussion can be found in Appendix A.1.

Both Coxeter and Freudenthal-Kuhn triangulations can be described as an arrangement of hyperplanes. They are related by an affine transformation. Let E be a finite set of vectors of properties of permutahedra are given in Appendix A. 1.1). We first recall that two complexes are dual if there is a bijection between their faces that inverses the inclusion relationships.

▶ Proposition 2. The star of a vertex in a CFK-triangulation is combinatorially dual to a permutahedron.

Since the facial structure of a permutahedron is fully described by ordered partitions, any simplex σ in a CFK-triangulation is characterized by a star that contains σ and by the ordered partition that specifies which simplex in the star is precisely σ. Since a simplex appears in several stars, we take the one that is centered at the lowest vertex of σ in the lexicographic order. This representation is called the permutahedral representation of a CFK-triangulation, see Figure 1. We further have:

PL-approximation of isomanifolds

We first recall sufficient conditions under which the PL-approximation M output by the algorithm faithfully reproduces the original isomanifold. These conditions are fully described in the full version of [START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF] and we simply state here the main results specialized to the case of CFK-triangulations.

We will say that f has bounded complexity if the three following quantities γ max , λ min and α max are positive and bounded.

γ max = max x∈T0 (max i |gradf i (x)|) λ min = min x∈T0 λ min (x), α max = max x∈T0 max i ∥Hes(f i )(x)∥ 2
where T 0 denotes the set of all σ ∈ T , such that (f i ) -1 (0) ∩ σ ̸ = ∅ for all i.

gradf i = (∂ j f i ) j denotes the gradient of component f i , for i ∈ [1, d -n],
Gram(∇f ) denotes the Gram matrix whose elements are ∇f i • ∇f j where • stands for the dot product. λ min (x) denotes the smallest absolute value of the eigenvalues of Gram(∇f (x)), 1 We can now restate the topological result of [START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF] :

Hes(f ) = (∂ k ∂ l f i ) k,
▶

Tracing isomanifolds

In this section, we describe an algorithm that computes a PL-approximation M of an isomanifold M. The algorithm has some similarity with the Marching Cube algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] but departs from it in two fundamental ways. First, because of the curse of dimensionality, we cannot afford to look at all the cells in the grid and need to limit the search to cells that are close to M. The problem of computing M can be naturally decomposed into two subproblems: locating the various components of M (i.e., finding at least one point in each connected component), and then tracing around each component, using the fact that the components are connected. This decomposition is used by various authors, see for example [START_REF] Wenger | Isosurfaces: geometry, topology, and algorithms[END_REF][START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF]. In this paper, we focus on the tracing problem, although we discuss very briefly 1 Because a Gram matrix is a symmetric square matrix, its eigenvalues are well defined and real. The second major difference with the original marching cube algorithm is to replace the usual cubical grid by a CFK-triangulation of the ambient space. Taking a CFK-triangulation instead of a grid is a major advantage in high dimensions that has been recognized in the pioneering works of Allgower, and Schmidt [START_REF] Allgower | An algorithm for piecewise-linear approximation of an implicitly defined manifold[END_REF] and of Dobkin et al. [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF], see also [START_REF] Min | Simplicial isosurfacing in arbitrary dimension and codimension[END_REF]. The novelty here is to use the data structure of Section 2.1 to represent a CFK-triangulation. As a consequence, we will keep two main advantages of using grids: very limited storage and fast basic operations.

Isomanifolds

Let f : R d → R d-n be a smooth (C 2 suffices) function, and suppose that 0 is a regular value of f , meaning that at every point x such that f (x) = 0, the Jacobian of f is non-degenerate.

Then the zero set of f is an n-dimensional manifold as a direct consequence of the implicit function theorem, see for example [START_REF] Duistermaat | Multidimensional Real Analysis I: Differentiation[END_REF]Section 3.5]. We further assume that f -1 (0) is compact.

As in [START_REF] Allgower | Estimates for piecewise linear approximations of implicitly defined manifolds[END_REF] we consider a triangulation T of R d . The function f is the linear interpolation of the values of f at the vertices if restricted to a single simplex σ ∈ T , i.e.

∀x ∈ σ : f (x) = v∈σ λ v (x)f (v), (1)
where the λ v are the barycentric coordinates of x with respect to the vertices v of σ. For any function g : R d → R d-n we write g i , with i = 1, . . . , d -n, for the components of g.

The PL-approximation is now defined as f -1 (0 and σ, it is an n-dimensional polytope denoted by C σ . The PL-approximation or mesh M of M is the polytopal cell complex obtained by gluing the polytopes C σ associated to all the simplices σ in T .

) = M. Locally, f | -1 σ (0) is

Manifold tracing algorithm

Let M be the zero set of some function f : R d → R d-n , and let M be the associated PL-approximation defined over a triangulation T of the ambient space R d . Both n, and d are known but arbitrary, and will be considered as parameters in the complexity analysis. We write k = d -n for the codimension of M. The algorithm will use for T a CFK-triangulation stored using the data structure from Section A.1. We assume that the manifold M, and the triangulation T satisfy the following genericity hypothesis:

▶ Hypothesis 8 (Genericity We note that this condition can be satisfied by an infinitesimal perturbation for isomanifolds.

This requires some explanation. We recall that the CFK-triangulation is a hyperplane arrangement, and up to translation there are a finite number of k-flats that contain all k-simplices in the CFK triangulation. Hypothesis 8 is not satisfied, if either the flat H σ is not linearly independent of these k-flats, or if H σ does intersect some (k -1)-flat in the CFK-triangulation. In the previous section, we have already seen that an infinitesimal perturbation ensures that H σ is n-dimensional. Because two affine spaces whose dimensions do not add up to the ambient dimension don't intersect with generically and two affine spaces whose dimensions add up to exactly the ambient dimension intersect in a single point, we see that genericity can be achieved by perturbing f infinitesimally. We further remark that, generically, any vertex of the PL-approximation M is the intersection point between a k-simplex σ of T with the n-flat H σ that interpolates f inside σ. The seed simplices are given as part of the input and we don't discuss in this paper the problem of their construction. We simply observe that they can be obtained by computing a critical point (e.g., a point with smallest x 1 -coordinate) on each connected component of the isomanifold, which reduces to finding a solution to a system of equations, on which a large body of literature exists. See for example [START_REF] Alexander | Solution of Equations and Systems of Equations: Pure and Applied Mathematics: A Series of Monographs and Textbooks[END_REF][START_REF] James | Iterative solution of nonlinear equations in several variables[END_REF][START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF] and also the discussion in Wenger's book [START_REF] Wenger | Isosurfaces: geometry, topology, and algorithms[END_REF]Section 8.4]. Once such a seed point has been computed, we simply translate and rotate the triangulation T so that the seed point coincides with the barycenter of a k-simplex of T and the intersection with the manifold is transversal as demanded by the genericity hypothesis (for numerical stability it is convenient if the angle between the tangent space of the manifold and the starting k-simplex is large, which is easy to ensure). If the distance between M and M is small enough, then M also intersects the same k simplex, see Section 2.2.

The algorithm is described as Algorithm 1. It takes as input the permutahedral representation of an ambient FKC-triangulation T and a seed k-simplex τ 0 of T . We assume that T satisfies the Genericity Hypothesis 8, which can be enforced by infinitesimal perturbations of f as discussed in Section 3.1.

The algorithm maintains the subset S of the simplices in T of dimension k that intersect M.

S is initialized with the seed simplex τ 0 and stored as a hash table so that we can decide in constant time if a given k-simplex belongs to S. Then, starting from τ 0 , we look at all its cofacets and consider all the facets of those cofacets that are not in S (i.e. they have not been considered yet). This can be done using a queue Q of candidate k-simplices. Each of these simplices is queried with the intersection oracle and, if it is found to intersect M, it is added to S if not already present. Upon termination, S contains all the k-dimensional simplices of T that intersect M. Each such intersection, which consists of a single point (by the Genericity hypothesis), is a vertex of M. Hence M0 is the vertex set of M.

Note that our algorithm essentially traverses the adjacency graph of the k and (k+1)-simplices of T that intersect M. It therefore identifies not only the set M0 of vertices of M, but also the edges joining two such vertices (associated to the cofacets of the k-simplices in S). By simply reporting those cofacets on the fly, the algorithm can output the 1-skeleton M1 of the n-dimensional polytopal cell complex M. The higher dimensional faces of M are the polytopes C τ = τ ∩ H τ for all the cofaces τ of the k-simplices of S. If needed, the full Hasse diagram of M can be computed from M0 . This can be done in an output sensitive manner by using the permutahedral representation of T and the algorithm of Section A.1 to compute cofaces by increasing dimensions.

Complexity analysis

We can easily bound the complexity of the manifold tracing algorithm as a function of the size of the output. We will now express the size of the output in terms of quantities that depend on the manifold, the ambient dimension d, and the resolution of the triangulation (the diameter D of a simplex) which bounds the density of the output sample, and the precision of the approximation. Our result holds for K-sparse manifolds, i.e. submanifolds whose intersection with any k-flat consists of at most K points. In practical situations, K is usually small and, in particular, K is a constant for algebraic isomanifolds of bounded degree.

N C ≤ K n! × d 2 √ d(d+2) 2 √ 2D n and N F K ≤ K n! × d 3 √ 2D n
where D is the diameter of a simplex of T .

We see that Coxeter triangulations lead to smaller samples than FK-triangulations by a factor of roughly 2 n . This will be confirmed experimentally (see Figure 3).

As noticed in Section 3.2, a simple variant of the algorithm can compute the full Hasse diagram of M in an output sensitive manner. The following lemma shows that the combinatorial complexity of M is of the same order as the combinatorial complexity as M0 .

▶ Proposition 11. The combinatorial complexity of M is |S| × ( 3 2 ) n (n + 1)!,

Dimensionality reduction

As seen from Proposition 10, the size |S| of the output of the algorithm, considered as a function of the resolution D of the triangulation, depends exponentially on n (which is to be expected), and only polynomially on d (which is fortunate). Nevertheless, the computing time of our algorithm and the size of the output depend on d. Removing the dependency on d in the time complexity is impossible since we need to evaluate a vector-valued function f at a number of points of R d , which takes Ω(d) time per evaluation. However, we will see that we can reduce the size of the mesh produced by our algorithm.

Examples of samples of M whose sizes depend on n but not on d, and lead to good Since its cardinality depends on d, it is not well separated and, in particular, not a D-net of M. If we are mostly interested in the output sample, we can easily sparsify it to obtain a D-net. However, by doing so, we will lose the combinatorial structure of the mesh.

We now show how to compute a D-dense sample of M of size independent of d, together with a mesh. Specifically, we will reduce dimensionality using a variant of the celebrated Johnson-Lindenstrauss lemma for manifolds. Doing so, we depart from our previous worstcase analysis by allowing some approximation factor ε and tolerate a guarantee that holds only with high probability.

▶ Theorem 13 (Johnson-Lindenstrauss lemma for manifolds [START_REF] Kenneth | Tighter bounds for random projections of manifolds[END_REF][START_REF] Verma | A note on random projections for preserving paths on a manifold[END_REF]). Pick any ε, η > 0,

and let d ′ = Ω n ε 2 log 1 ε + 1 ε 2 log Γ η ,
where Γ is a quantity that depends only on intrinsic properties of M. Let Φ be the projection on a random affine subspace of dimension d ′ . Then,

with probability > 1 -η, for all x, y ∈ M, we have (1 -ε) d ′ d ≤ ∥Φx-Φy∥ ∥x-y∥ ≤ (1 + ε) d ′ d . Let Ψ = d d ′ Φ. By the theorem, the image Ψ(M) of M is a submanifold of dimension n embedded in R d ′ .
One can now run the manifold tracing algorithm in R d ′ to sample, and mesh Ψ(M). The algorithm works as described before except that we need another oracle that, given a (d

′ -n)-simplex σ of the CFK-triangulation of R d ′ , decides whether its inverse image Ψ -1 (σ) intersects M or not. Note that Ψ -1 (σ) is a (d -d ′ )-dimensional flat strip
(that is the product of a face and an affine subspace) in R d , and that the complexity of this new oracle is the same as the complexity of the basic intersection oracle, i.e. polynomial in d.

Due to the scaling factor d/d ′ , the resolution of the triangulation in the low dimensional space R d ′ has to be scaled by the same factor if one wants to satisfy a given sampling density on M. Since the geometry of the manifold is also scaled in the same way [START_REF] Eftekhari | What happens to a manifold under a bi-Lipschitz map?[END_REF], the analysis of the algorithm will be unchanged. Proposition 10 then shows that the size of the output sample does not depend on d but only on n and D for fixed ε, and η. Moreover, since the complexities of the projection and of the new oracle are polynomial in d, Proposition 9

implies that the overall complexity is still polynomial in d.

Isomanifolds with boundary, and isostratifolds

The case of isomanifolds with boundary and, more generally, of isostratifolds can be handled in very much the same way. By an isomanifold of dimension n with boundary, we mean that, on top of a function f : R d → R d-n , we are given another function f ∂ : R d → R, and the set

we consider is M = f -1 (0) ∩ f -1 ∂ ([0, ∞)). We note that ∂M = f -1 (0) ∩ f -1 ∂ (0).
Similarly to (1), we also define f∂

| τ (x) = v∈σ λ v (x) f ∂ (v).
We write f for the (global)

piecewise linear function that coincides with f | τ on each τ of T , and f∂ for the (global)

piecewise linear function that coincides with f∂ | τ on each τ of T . We note that the piecewise linear approximation of the boundary

∂ M = f -1 ∂ (0) ∩ f -1 (0) is a subset of f -1 (0), i.e. the
piecewise linear approximation of the manifold ignoring the boundary. The piecewise linear approximation M of the manifold with boundary consists of the following cells:

For each τ of T , such that f∂ | τ is positive on τ , and ( f

| τ ) -1 (0) ∩ τ ̸ = ∅, we add ( f | τ ) -1 (0) ∩ τ . For each τ of T , such that ( f | τ ) -1 (0) ∩ τ ̸ = ∅, and ( f∂ | τ ) -1 (0) ∩ τ ̸ = ∅, we add ( f | τ ) -1 (0) ∩ ( f∂ | τ ) -1 ([0, ∞)) ∩ τ .
We will assume that the Genericity Hypothesis 8 holds for both M, and ∂ M.

We can now adapt the algorithm of Section 3.2 as follows. In addition to reporting the set S k of k-faces of the triangulation T that intersect M, the algorithm will also report the set The three rules above together with the permutahedral representation of T provide a way to construct the 1-skeleton of M on the fly. The total cost is output sensitive. If needed, the entire combinatorial structure of M can be computed by traversing the full triangulation T .

The above construction generalizes easily to arbitrary isostratifolds. Isostratifolds are stratified spaces that are defined by equations and inequalities. An example of such a stratifold is an octant of the sphere in R 3 that can be defined by as

x 2 + y 2 + z 2 -1 = 0,
x ≥ 0, y ≥ 0, and z ≥ 0. We compute the 1-skeleton of M and construct a graph whose nodes are the simplices of dimensions k, k + 1, ..., d that intersect the strata of dimension n, n -1, ..., 0.

Experimental results

The data structure of Section A.1 and the algorithm of Section 3 have been implemented in C ++ . The code is robust and fast and will be released in the GUDHI library [START_REF]GUDHI Project[END_REF]. Full detail on the implementation, including the implementation of the oracle, can be found in [START_REF] Kachanovich | Meshing submanifolds using Coxeter triangulations. Theses, COMUE Université Côte d'Azur[END_REF].

In this section, we explore the dependency of our C ++ implementation of the data structure for the ambient CFK-triangulation, and of the manifold tracing algorithm on the properties of the triangulation, and of the input manifold.

Performance of the algorithm

We show the performance of our implementation of the manifold tracing algorithm for various ambient and intrinsic dimensions in Figure 2. In Figure 3, we can see that using

Coxeter triangulation is beneficial in practice as it produces a smaller output in less time (see Proposition 10).

In Figure 4, we present a PL approximation of a two-dimensional Clifford torus without boundary embedded in R 10 built by the manifold tracing algorithm. The torus has been rotated and translated in R 10 so that the coordinate axes do not play any special role. Note that there is no C 2 isometric embedding of the Clifford torus in R

(x 2 1 + x 2 2 + x 2 3 -0.8) 2 -0.4 (x3 -1) 2 -2x 2 1 (x3 + 1) 2 -2x 2 2 =

Manifolds with boundary

The algorithm has been adapted to handle submanifolds with boundary and surfaces with a piecewise smooth boundary, see Section 3.5. In Figure 5, we present the mesh obtained by our algorithm on a portion of a flat torus embedded in R 4 , and cut by a hypersphere. The torus has been rotated and translated in R 4 so that the coordinate axes do not play any special role.

An application in algebraic geometry

We also applied our algorithm to a more complicated example of interest in algebraic geometry [START_REF] Alvarez | Dynamique et topologie du feuilletage de Jouanolou[END_REF] where an active field of research is to understand the geometry and topology of various projective varieties. Projective varieties are isomanifolds defined by polynomial equations in the complex projective space CP d = (C d+1 \ 0)/C * of complex dimension d.

One such example is the complex one-dimensional curve (that is a real dimensional surface)

given by the equation z 2 1 z2 + z 2 2 z3 + z 2 3 z1 = 0 in CP 2 , where z denotes the conjugate of the complex number z.

To be able to apply our algorithm, we first need to pass from homogenous coordinates [z 1 : . . . :

z d+1 ] on CP d to affine coordinates [z ′ 1 : • • • : z ′ i-1 : 1 : z ′ i+1 : • • • : z ′ d+1
] by picking the ith coordinate to be equal to 1, that is z ′ j = z j /z i . Given some homogenous coordinates [z 1 : . . . : z d+1 ], we can choose the ith coordinate to be set to 1 to be the coordinate whose absolute value is the largest, so that CP d can be written as the union of the d + 1 sets

[z ′ 1 : • • • : z ′ i-1 : 1 : z ′ i+1 : • • • : z ′ d+1 ] | |z ′ j | ≤ 1
, with the boundaries of these sets identified.

Writing z ′ j = x j + iy j these sets are (seen as real sets) identical to the domain of R 2d

D i = {(x 1 , y 1 , . . . , x i-1 , y i-1 , x i+1 , y i+1 , . . . , x d+1 , y d+1 ) | x 2 j + y 2 j ≤ 1}.
Let f be a homogenous polynomial in d + 1 complex variables and their complex conjugates.

For each i, we can fix the ith coordinate to be 1. Writing each variable in terms of its real This principle generalizes to varieties of higher codimension, that is to varieties defined by a number of homogenous polynomials f 1 , . . . , f d-m .

Figure 6 The three triangulated surfaces as discussed in the example of z 2 1 z2 + z 2 2 z3 + z 2 3 z1 = 0 in CP 2 after projection from R 4 to R 3 .

We illustrate the above construction on the above equation

z 2 1 z2 + z 2 2 z3 + z 2 3 z1 = 0 in CP 2 .
By passing to affine coordinates, we recover z 2 1 z2 + z 2 2 + z1 = 0, z 2 1 + z3 + z 2 3 z1 = 0, and z2 +z 2 2 z3 +z 2 3 = 0. By expanding z 1 = x 1 +iy 1 , z 2 = x 2 +iy 2 , and z 3 = x 3 +iy 3 , we find two real equations for each of the complex equations. We give those corresponding to z 2 1 z2 +z 2 2 + z1 = 0, the other equations being symmetric. For this complex equation, we get the real equations

x 1 + x 2 1 x 2 + x 2 2 -x 2 y 2 1 + 2x 1 y 1 y 2 -y 2 2 = 0 and -y 1 + 2x 1 x 2 y 1 -x 2 1 y 2 + 2x 2 y 2 + y 2 1 y 2 = 0 in R 4 .
The domain D 3 is in this case determined by the equations

x 2 1 + y 2 1 ≤ 1 and x 2 2 + y 2 2 ≤ 1.
Hence we find a surface in R 4 with a piecewise smooth boundary. The result provided by our algorithm is shown in Figure 6. For visualization purposes, we show the three surfaces separately and projected from R 4 to R 3 .

S o C G 2 0 2 1 Conclusion and open questions

We have presented an efficient, practical and provably correct algorithm to compute the PL-approximation of an isomanifold of any dimension and codimension. Since isomanifolds are a special type of manifolds, it is tempting to see if our algorithm extends to general smooth submanifolds of R d .

The manifold tracing algorithm itself is quite general and works for any submanifold as soon as we provide a seed point and an oracle that can determine whether a k-simplex of the ambient triangulation intersects M or not. In this general setting, the simple algorithm described above is sufficient to compute a PL-approximation of the manifold and satisfies the bounds given in Section 3.

However, this is not enough to obtain guarantees on the geometric and topological quality of the output mesh. Such guarantees can be obtained by slightly perturbing the ambient Coxeter triangulation of type Ãd so that the following conditions are satisfied:

1. All k-dimensional faces τ in T , with k ≤ d -n -1, are far enough from M.

2.

The longest edge length of T is upper bounded and its smallest thickness is lower bounded.

Under these conditions, Algorithm 1 will output a PL-approximation that is topologically equivalent and close in Hausdorff distance to the input manifold [START_REF] Boissonnat | Triangulating submanifolds: An elementary and quantified version of whitney's method[END_REF]. However, the perturbation scheme of [START_REF] Boissonnat | Triangulating submanifolds: An elementary and quantified version of whitney's method[END_REF] perturbs (in the worst case) all the simplices of T of dimension less than the codimension d -n that are incident on a vertex (in a neighbourhood of M).

Since there are exponentially many such simplices, such methods have a complexity that depends exponentially on the ambient dimension d, and have not proved useful in practice except in some simple cases. It remains open whether general smooth manifolds embedded in R d can be triangulated in time polynomial in d as we were able to do here in the special case of isomanifolds.

A A compact data structure for high-dimensional

Coxeter-Freudenthal-Kuhn triangulations

Subdivisions of Euclidean space are a major tool to efficiently answer geometric queries, compute approximation of shapes or solve optimization problems. Among the most widely used subdivision schemes are grids and triangulations. Both are subject to the curse of dimensionality and their combinatorial complexity depends exponentially on the dimension of the space. Triangulations are most flexible since their vertex set can be any set of points.

Differently, uniform grids depend only on the space but not on a given data set. The rigidity of the grid structure has a major advantage: the grid, although of exponential size, need not be represented explicitly and basic operations like locating a point or computing faces or cofaces of a given cell in the grid can be done without storing an explicit representation of the grid. In fact, the representation can be entirely implicit. This is clearly impossible with general triangulations with arbitrary vertex sets.

The question of designing efficient data structure for triangulations and more general simplicial complexes led to interesting developments recently. On one hand, one can take advantage of the fact that special types of simplicial complexes allow compact representations. Most notably, flag complexes (including the celebrated Vietoris-Rips complex) can be represented by their 1-skeleton (or graph) and higher dimensional faces can be retrieved by computing the cliques of the graph. One can also represent a simplicial complex by its blockers, i.e. the simplices that do not belong to the complex but whose facets do [START_REF] Attali | Efficient data structure for representing and simplifying simplicial complexes in high dimensions[END_REF].

On a different front, data structures have been proposed to efficiently store general simplicial complexes such as the simplex tree [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF] that uses a trie to store the faces of all dimensions, or the Simplex Array List [START_REF] Boissonnat | Building efficient and compact data structures for simplicial complexes[END_REF] that represents only the maximal faces, which allows an exponential saving in storage since a simplex has exponential complexity. Nevertheless, due to their generality and the fact that the represented complexes don't have any prespecified symmetry, the data structures cannot compete with grids in terms of size and efficiency.

The Coxeter-Freudenthal-Kuhn triangulations in this paper form a middle ground, they form a special class of triangulations of R d that have a high regularity. The data structure for such triangulations that can be as compact as for grids.

The Coxeter-Freudenthal-Kuhn triangulations we consider combine two classes of triangulations with different origins and names. The two foundational works are due to Coxeter [START_REF] Harold | Discrete groups generated by reflections[END_REF] and Freudenthal [START_REF] Freudenthal | Simplizialzerlegungen von beschrankter flachheit[END_REF]. Coxeter triangulations derive from geometric group theory, in particular affine Weyl groups, while Freudenthal triangulations (also called Kuhn triangulations) are combinatorial in nature. Nevertheless, both triangulations are the same up to a linear transformation, as remarked in [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF] and fully proved in this paper. This allows us to combine the nice geometric properties of Coxeter triangulations of type Ãd with the simple combinatorial definitions of the Freudenthal-Kuhn triangulation, and its connection to permutahedra.

Coxeter triangulations of type Ãd are geometrically attractive because each simplex is very well shaped (large volume compared to longest edge length), and all d-simplices are identical up to reflections.

Although these triangulations do not depend on a given data set, they proved to be very useful in to interpolate multivariate multivalued functions or to mesh geometric shapes embedded in high dimensional spaces. Freudenthal-Kuhn triangulations have been known in Applied Mathematics [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF][START_REF] Eaves | A course in triangulations for solving equations with deformations[END_REF][START_REF] Michael | The computation of fixed points and applications[END_REF], and Coxeter triangulations have been used by Dobkin et al. [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF] to trace curves in high dimensions and are good candidates to trace manifolds of any codimension [START_REF] Min | Simplicial isosurfacing in arbitrary dimension and codimension[END_REF]. They have also been used in the context of Topological Data Analysis [START_REF] Choudhary | Polynomial-sized topological approximations using the permutahedron[END_REF].

In Appendix A.1, we study these triangulations. This section recalls and extends to arbitrary dimensions several results that were disseminated in many different places which are sometimes difficult to access and in different languages (see among others [START_REF] Freudenthal | Simplizialzerlegungen von beschrankter flachheit[END_REF][START_REF] Michael | The computation of fixed points and applications[END_REF][START_REF] Eaves | A course in triangulations for solving equations with deformations[END_REF][START_REF] Maes | On the permutahedron and the quadratic placement problem[END_REF][START_REF] Ziegler | Lectures on Polytopes[END_REF][START_REF] Rennie | On Stirling numbers of the second kind[END_REF][START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF]).

Based on these results, we introduce in Appendix A.2 a very compact data structure that implicitly stores the full facial structure of such triangulations. The data structure allows to locate a point in the triangulation and to retrieve the faces or the cofaces of a simplex of any dimension in an output sensitive way.

The Data Structure has been implemented and fully tested. Section 4 reports on experimental results and demonstrates that the data structure is highly practical and can handle triangulations of high dimensional spaces. Using our data structure will allow to extend the applicability of the methods based on such triangulations and to significantly improve their performance. It appears to be especially useful to trace low dimensional manifolds embedded in high dimensional spaces as encountered in statistics, dynamical systems, econometrics, or mechanics [START_REF] Chen | Solution manifold and its statistical applications[END_REF][START_REF] Michael | The computation of fixed points and applications[END_REF][START_REF] Newman | A survey of the marching cubes algorithm[END_REF].

A.1 Coxeter-Freudhental-Kuhn triangulations

Freudenthal-Kuhn triangulations are combinatorial structures that come from a specific triangulation of the d-cube. Their connections to permutahedra is at the heart of our data structure. Coxeter triangulations, we introduce in Section A.1.3, have a different flavour and come with very nice geometric properties. Since both types of triangulations are the same up to an affine transformation, as first noted by Dobkin et al [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF], they have the same combinatorial structure and our data structure will be able to handle both of them.

Although most ideas in this section were known previously, we give full proofs of the results that were not explicitly mentioned or not proved in full generality in the literature.

A.1.1 Permutahedra

We write [i] = {1, . . . , i}, and [i, j] = {i, . . . , j}.

▶ Definition 14 (Permutahedron). A d-permutahedron is a d-dimensional polytope, which is the convex hull P of all points in R d+1 , the coordinates of which are permutations of

[d + 1].
Formally, this convex hull can be written as:

P = conv{(σ(1), . . . , σ(d + 1)) ∈ R d+1 | σ ∈ S d+1 },
where S d+1 denotes the set of permutations of [d + 1]. If there is a need to distinguish permutahedra of different dimension we write P(n) the permutahedron of dimension n -1.

P is at most d-dimensional since all its vertices lie on the hyperplane of equation

d+1 i=1 x i = d(d + 1) 2 .
Moreover, it can be shown that there are d + 1 affinely independent vertices in P, proving that P is exactly d-dimensional (see for example [START_REF] Maes | On the permutahedron and the quadratic placement problem[END_REF]Lemma 3.4]). The facial structure of P is best described in terms of ordered partitions [START_REF] Ziegler | Lectures on Polytopes[END_REF]. Refer to Figure 7.

2 2(d+1) log(i+1) .
The following three corollaries seem to be new.

▶ Corollary 19. The number p 0,i of vertices of a i-face of a d-permutahedron is at most (i + 1)!, and at least 2 min(i,d-i) .

The proof of Corollary 19 is based on:

▶ Lemma 20 (Lemma 3.11 of [START_REF] Maes | On the permutahedron and the quadratic placement problem[END_REF]). The face of a permutahedron corresponding to an ordered partition ω = (ω 1 , . . . , ω l+1 ) is combinatorially Using the inverse argument, the lower bound is obtained when each part in the ordered partition are as small as possible that is all parts have almost equal size. In this case,

P(|ω 1 |) × • • • × P(|ω l+1 |),
|ω p | ≥ ⌊ d+1 l+1 ⌋, so that l+1 p=1 (|ω p |!) ≥ d + 1 l + 1 ! l+1
More accurately, let r ′ be the remainder of d + 1 after division by l + 1, that is r ′ = d + 1 mod l + 1 then:

l+1 p=1 (|ω p |!) ≥ d + 1 l + 1 ! l-r ′ +1 d + 1 l + 1 + 1 ! r ′
We now distinguish two cases: We have

If ⌊ d+1 l+1 ⌋ ≥ 2,
p i,j ≤ 1 2 min(i,d-i)
c v = j i |F 0,σ | ≤ j i (j + 1)!.
In addition, by Corollary 19, we have

τ ∈Fi,σ p τ ≥ 2 min(i,d-i) |F i,σ |
The inequality follows since σ is any j-face of the d-permutahedron. ◀

A.1.2 Freudenthal-Kuhn triangulation

The Freudenthal-Kuhn (FK for short) triangulation is obtained from the d-grid, i.e. the unit cubical tessellation of R d that consists of copies of the unit d-cube along the integer lattice Z d . By triangulating each d-cube in the grid in an appropriate way to be described now, we obtain the FK-triangulation of R d . The results and definitions below were known to Freudenthal [START_REF] Freudenthal | Simplizialzerlegungen von beschrankter flachheit[END_REF], Todd [START_REF] Michael | The computation of fixed points and applications[END_REF], or Eaves [START_REF] Eaves | A course in triangulations for solving equations with deformations[END_REF], mainly for top dimensional simplices and in different guises. We combine these results and extend to simplices of arbitrary (co)dimension, where necessary.

▶ Definition 23. Let x ∈ R d , and write z i = x i -⌊x i ⌋. We denote by e 1 , . . . , e d the basis vectors and introduce, for reasons that will be clear later, the extra vector

e d+1 = - d i=1 e i .
We introduce the convention that z d+1 = 0. We associate to x the ordered partition ω = (ω 1 , . . . , ω l+1 ) of [d + 1] where the ω i are obtained by sorting the z i in decreasing order.

Specifically, with ω i = {ω i (1), . . . , ω i (m i )}, we have (see Figure 8):

1 >z ω1(1) = • • • = z ω1(m1) > • • • > z ω l (1) = • • • = z ω l (m l ) > z ω l+1 (1) = • • • = z ω l+1 (m l+1 ) = 0. (2) {1, 2, 3} {1}, {2}, {3} {2}, {1}, {3} { 1 , 2 } , { 3 } {2}, {1, 3} {1}, {2, 3}
v 0 = (⌊x 1 ⌋, . . . , ⌊x d ⌋), v i = v i-1 + j∈ωi e j i = 1, . . . , l. ( 3
)
Then x is a point in the relative interior of σ if and only if

z i = x i -⌊x i ⌋, i = 1, . . . , d + 1
(with, as above, z d+1 = 0), satisfy (2).

Proof. Because the whole problem is translation invariant, we assume that v 0 = 0 without loss of generality, so that the expressions are shorter. Using barycentric coordinates, z ∈ σ can be written as

z = l i=0 λ i v i = l i=0 λ i i k=1 j∈ωi e j =λ l k∈ω l e k + (λ l + λ l-1 )   k∈ω l-1 e k   + • • • + (λ l + • • • + λ 1 ) k∈ω1 e k . ( 4
)
Here the

λ i > 0, i ∈ [0, l],
l i=0 λ i = 1, are the barycentric coordinates of z in σ. We have

α ω l (1) = • • • = α ω l (m l ) = λ l . . . α ω1(1) = • • • = α ω1(m1) = λ l + • • • + λ 1 (5)
From (4), we see that α ωi(j) is the ω i (j)th coordinate of z, denoted by z ωi(j) , while all coordinates z ω l+1 (1) , . . . , z ω l+1 (m l+1 ) are zero, because e ω l+1 (i) does not occur in (4), for all i.

Moreover, because λ l + • • • + λ i > λ l + • • • + λ i-1 , we see that (2) is satisfied.
Conversely, given a point z such that its coordinates satisfy (2), we can read of its barycentric coordinates with respect to the v i , as defined by (3), from [START_REF] Armstrong | Groups and Symmetry[END_REF]. ◀ ▶ Theorem 25. The equivalence classes of the points of R d with a same ordered partition are the simplices of a triangulation of R d called the FK-triangulation (see Figure 8).

Proof. Lemma 24 implies that:

Any face of a simplex in the FK-triangulation also lies in the FK-triangulation.

The intersection of two simplices in the FK-triangulation also lie in the FK-triangulation.

For any point x ∈ R d , there is a unique simplex σ such that x lies in the relative interior of σ. Indeed, x has uniquely defined barycentric coordinates with respect to the vertices of σ, and thus is mapped to a unique point in σ.

Hence the partition we have defined is a well-defined triangulation of R d . ◀ ▶ Remark 26. We note that, by construction, v 0 in Lemma 24 is the smallest vertex of σ in the lexicographical order. Lemma 24 also implies an observation of Freudenthal [START_REF] Freudenthal | Simplizialzerlegungen von beschrankter flachheit[END_REF]: all d-simplices in the FK-triangulation can be described by monotone paths along the edges of the cube from vertex (0, . . . , 0) + v 0 to vertex (1, . . . , 1) + v 0 . Conversely, any monotone path along the edges of the cubes from (0, . . . , 0) + v 0 to (1, . . . , 1) + v 0 gives a simplex in the FK-triangulation.

A.1.3 CFK-triangulations

Freudenthal-Kuhn triangulations are closely related to Coxeter triangulations of type Ãd [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] and both are arrangements of hyperplanes as demonstrated now.

Let E be a finite set of vectors of R d , and consider the set of hyperplanes

H E = {x ∈ R d | ⟨x, u⟩ = k, u ∈ E, k ∈ Z}.
In generalization of the theory of Coxeter triangulations, we call the set E the set of roots for historic reasons, as mentioned below.

These hyperplanes partition R d in a cell complex called the arrangement of the hyperplanes.

We denote it by H E .

▶ Lemma 27. The Freudenthal-Kuhn triangulation is the hyperplane arrangement

H E F K associated to the set of vectors E F K = {e 1 , . . . , e d } ∪ {u i,j = e j -e i | 1 ≤ i < j ≤ d}.
Proof. Thanks to Lemma 24, we know that if x ∈ σ, where σ is a simplex of dimension less than d, there is at least one equality in (2) on top of z d+1 = 0. That is x iṽi 0 = x jṽj 0 or

x iṽi 0 = 0 for some i, j ̸ = d + 1. Note that ṽi 0 , ṽj 0 ∈ Z. The converse direction of Lemma 24

gives that if x iṽi 0 = x jṽj 0 or x iṽi 0 = 0 for some i, j ̸ = d + 1 for x ∈ R d , then there is a simplex σ of dimension strictly less than d in the FK-triangulation such that x ∈ σ. ◀

Observe that the norms of the vectors in E F K are either 1 or √ 2. By definition, this implies that the distance between the two hyperplanes ⟨x, u⟩ = k, and ⟨x, u⟩

= k + 1, where u ∈ E F K , is either 1 or 1/ √ 2.
Let H be the hyperplane of R d+1 of equation ⟨x, 1⟩ = 0 where 1 is the vector of R d+1 whose coordinates are all 1. We now define a linear map µ from R d to H by showing how it acts on the basis of R d : µ(e i ) = r 1,i = i j=1 s j , where s i = e i -e i+1 , i = 1, ..., d. The vectors s j are called simple roots and play an important role in algebra. We refer to [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] for more information.

▶ Lemma 28. µ maps E F K bijectively onto the set E C defined as

E C = r i,j = j l=i s l = e i -e j+1 1 ≤ i ≤ j ≤ d
Proof. The vector µ(e i ) = r 1,i lies in E C , by definition. For u i,j ∈ E F K , with i < j, we see that

µ(u i,j ) = µ(e j -e i ) = µ(e j ) -µ(e i ) = r 1,j -r 1,i = j l=1 s l - i l=1 s l = j l=i+1 s l = r i+1,j .
Hence µ(u i,j ) lies in E C . By reading the previous calculation backwards, we see that µ -1

maps each r ∈ E C to a vector in E F K . ◀
Observe that all vectors in E C have length √ 2. By definition, this implies that the distance between the two hyperplanes ⟨x, u⟩ = k, and ⟨x, u⟩

= k + 1, where u ∈ E C , is 1/ √ 2.
The image by µ of the Freudenthal-Kuhn triangulation is a triangulation of R d which is the arrangement H E C associated to the set of vectors E C . This triangulation is called the Coxeter triangulation of type Ãd of R d . By definition, it has the same combinatorial structure as the FK-triangulation. In addition, it has remarkable geometric properties [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF][START_REF] Choudhary | Coxeter triangulations have good quality[END_REF]. First, it is a non-degenerate Delaunay triangulation, and its dual complex is a Voronoi diagram.

Moreover, its simplices have an exceptionally large thickness (the ratio of the smallest altitude of a simplex over its diameter or longest edge length).

A.2 Data Structure

We introduce our data structure in this section. We first consider the case of FK-triangulations in Sections A.2.1 and A.2.2. The extension to CFK-triangulations in Section A.3 is straightforward, since all these triangulations have the same combinatorial structure.

A.2.1 Permutahedral representation of FK-triangulations

Cycles and the permutahedron. In Remark 26 we have seen that simplices can be described by monotone paths (increasing coordinates) along the edges of the cube. As observed by Eaves [START_REF] Eaves | A course in triangulations for solving equations with deformations[END_REF], these monotone paths can be made into a cycle using the extra vector e d+1 = -

d i=1 e i
because by construction

v 0 = v l + i∈ω l+1 e i ,
with ω as in Definition 23. Because it is a cycle, we can take any vertex of the cycle as a starting point, which means that v 0 no longer has a special role as a starting point of a monotone edge walk. A cycle can now be represented by an ordered partition of [d + 1], for which it is not longer necessary that d + 1 lies in ω l+1 , and an (arbitrary) starting point. We now formalize these general cyclical paths:

▶ Definition 29 (Permutahedral representation). Let (v 0 , ω) ∈ Z d × OP l+1 [d].
To this pair we associate a simplex

σ ω = {v 0 = v ω 0 , v ω 1 , . . . , v ω l } with v ω i = v ω i-1 + i∈ωi e i i = 1, . . . , l. ( 6
)
We say that (v 0 , ω) is the permutahedral representation of the simplex σ ω . If d + 1 ∈ ω l+1 , we say that (v 0 , ω) is the canonical permutahedral representation of σ ω . In this case, σ ω is a simplex in the FK-triangulation in the cube of which v 0 is the minimal vertex with respect to the lexicographical order, as we have seen above.

J.-D. Boissonnat, S. Kachanovich, and M. Wintraecken

XX:29

In Lemma 32, and Proposition 33, we will see that, more generally,

{(v 0 , ω) | ω ∈ OP [d + 1]}
is the star of v 0 in the FK-triangulation, where we identify simplices with their permutahedral representations.

▶ Definition 30 (Cyclic shifts). Let (v 0 , ω) be a permutahedral representation. We define the cyclic shift of (v 0 , ω) of length k to the left as (v ′ 0 , ω ′ ), where

v ′ 0 = v 0 + k j=1 i∈ωj e i ω ′ j = ω (j+k-1) mod (l+1)+1 . ( 7
)
Here we use the convention that the sum from 1 to 0 is empty. We write

(v ′ 0 , ω ′ ) = (v 0 , ω) ⊕ k.
▶ Lemma 31. The cyclic shift (v ′ 0 , ω ′ ) = (v 0 , ω) ⊕ k defines the same simplex as (v 0 , ω).

Proof. Follows by inserting [START_REF] Tom | Equilibrium triangulations of the complex projective plane[END_REF] in [START_REF] Attali | Efficient data structure for representing and simplifying simplicial complexes in high dimensions[END_REF]. ◀

We now prove that the all permutahedral representations for a fixed v 0 , form the star of v 0 . This is a crucial property that will be used to efficiently compute faces, and cofaces, and traverse the triangulation.

▶ Lemma 32. The set {(v 0 , ω) | ω ∈ OP [d + 1]}, where OP [d + 1] is the set of all ordered partitions of [d + 1],
gives all the simplices in the star of v 0 in FK-triangulation.

Proof. Let (v 0 , ω), with ω ∈ OP l+1 [d + 1], be such that d + 1 ∈ ω k . Let (v ′ 0 , ω ′ ) = (v 0 , ω) ⊕ (l -k + 1)
. By Definition 30, and Lemma 31, (v 0 , ω), and (v ′ 0 , ω ′ ) represent the same

simplex. Moreover d + 1 ∈ ω ′ l+1 , that is (v ′ 0 , ω ′ ) is a canonical permutahedral representation.
This implies that (v ′ 0 , ω ′ ) lies in the FK-triangulation by Lemma 24, and Theorem 25.

Conversely, suppose that (v ′ 0 , ω ′ ) is the canonical permutahedral representation of a simplex in the star of v 0 , that is there is some k such that v ′ k = v 0 , with v ′ k as in [START_REF] Allgower | An algorithm for piecewise-linear approximation of an implicitly defined manifold[END_REF]. Then (v 0 , ω) = (v ′ 0 , ω ′ ) ⊕ k is also a permutahedral representation of the same simplex. ◀ Faces. From (6) it is clear that merging two consecutive parts in the ordered partition ω = (ω 1 , . . . , ω l+1 ) corresponds to removing a vertex from the simplex, that is taking a facet.

Here we stress that we allow to merge ω 1 , and ω l+1 , but in that case we have to change the base point of the cycle to v 0 + l∈ω1 e l to obtain the canonical representation. For example, when looking at the two dimensional example in Figure 1, we see that the edges that contain y in the red triangle with permutahedral representation (y, ({1}, {2}, {3})) are (y, ({1, 2}, {3})), and (y, ({1}, {2, 3})). The third edge of the red triangle is (y ′ , ({2}, {1, 3})).

Generally, given an ordered partition ω in l + 1 parts all (l -j)-faces can be found by merging j consecutive parts in ω (for example merging ω 1 with ω 2 , and ω 3 with ω 4 ), where we allow ω l+1 to merge with ω 1 , but in this case we again need to change the base point to obtain the canonical representation.

Because the combinatorial structure of the faces is compatible with the permutahedron, Lemma 32 immediately gives the following result. We recall that two complexes are dual if there is a bijection between the faces that inverses the inclusion relationships (see for example This can be found in [21, Chapter 21, Section 3.F], see also [START_REF] Choudhary | Polynomial-sized topological approximations using the permutahedron[END_REF]. We note that Coxeter triangulations of type Ãd are combinatorially equivalent to FK-triangulations as discussed below in Section A.1.3. In the appendix we give a new proof that is more direct than the one in [21].

Duality. We can associate to a FK-triangulation T its dual complex T * . Since T is a simplicial complex, T * is a simple complex, that is a cell complex whose faces are all simple polytopes [START_REF] Boissonnat | Algorithmic Geometry[END_REF]. By Proposition 33, each d-dimensional cell of T * is a d-permutahedron. Face computation. Let σ be an l-simplex whose canonical permutahedral representation is (v 0 , ω), where ω is an ordered partition of [d + 1] into l + 1 parts. The computation of all k-faces of σ goes as follows. We use Ehrlich's subset generation algorithm [START_REF] Ehrlich | Loopless algorithms for generating permutations, combinations, and other combinatorial configurations[END_REF] to compute all the subsets of k + 1 elements from {v 0 , . . . , v l }. Let τ = {v m0 , . . . , v m k } be such a subset.

A.2.2 Basic operations on FK-triangulations

τ is a k-face of σ. We then compute the canonical permutahedral representation of all those k-faces τ .

We first sort the m i so that m 0 < • • • < m k using counting sort. Then, the canonical permutahedral representation (ṽ ′ 0 , ω ′ ) of τ is found by merging consecutive parts of ω so as to obtain k + 1 parts as follows :

v ′ 0 = v m0 = v 0 + j∈ω1 e j + . . . + j∈ωm 0 -1 e j ω ′ i = ω mi-1 ∪ . . . ∪ ω mi-1 for i ∈ {1, . . . , k} ω ′ k+1 = (ω 1 ∪ . . . ∪ ω m0-1 ) ∪ (ω m k ∪ . . . ∪ ω l+1 ).
The complexity of computing all subsets of k + 1 vertices of σ using Ehrlich's algorithm Coface computation. Computing the faces of a simplex σ consists in coarsifying its ordered partition. The computation of cofaces is the reverse. Here we refine the ordered partition.

Specifically, if σ is a k-simplex represented by its canonical permutahedral representation (v 0 , ω), and we want to compute its l-cofaces, we need to compute all refinements of ω into l + 1 parts.

More precisely, we need to subdivide each ω i in a i ≤ |ω i | subparts so that

k+1 i=1 a i = l + 1.
This can be done in time proportional to the number k + 1 of the generated subparts. We then need to consider all the permutations of these subparts since we are interested in ordered partitions. Using known algorithms by Walsh [START_REF] Walsh | Loop-free sequencing of bounded integer compositions[END_REF], and Ruskey and Savage [START_REF] Ruskey | Gray codes for set partitions and restricted growth tails[END_REF], we can compute all the ordered partitions associated to the l-cofaces of σ in time proportional to the number of such cofaces. We thus obtain all the permutahedral representations (v 0 , ω ′ ) of all the l-cofaces of σ.

It is important to notice that all cofaces of σ have v 0 as a vertex. However v 0 is not necessarily the minimal vertex of some of the computed cofaces. We thus have to identify the minimal vertex of each computed coface, and use cyclic shifts (as in Lemma 32)

s = p d-l,d-k ≤ 1 2 min(l,d-l) d -k d -l (d -k + 1)!
is the number of l-cofaces of a k-simplex in the FK-triangulation. The space complexity of the algorithm is O(d).

A.3 Data structure for CFK-triangulations

We store a CFK-triangulation as follows. The combinatorial structure of the triangulation is given through the canonical permutahedral representation of its simplices, and the algorithms from Section A.

A.4 Experimental results for the datastructure

The data structure and the basic operations have been implemented in C++ and are currently under review to be integrated in the GUDHI library. We report on the execution time of the face and coface generation algorithms for the FK-triangulations.

In Tables 1234, we consider an ambient space of moderate dimension d = 30 and compute the higher dimensional faces of various high dimensional simplices, of dimensions ranging from 22 to 30.

Each entry in Table 1 corresponds to the total time in milliseconds of computing all the k-dimensional faces of a set of l-dimensional simplices in R 30 . The l-dimensional simplices are picked at random in the triangulation and the results are averaged over 1 000 simplices.

Note that the time 11 904.7ms is the time of computing all 5 852 925 faces of dimension 22 of a simplex of dimension 30. In Tables 3 and4, we present analogous tables for the coface computation algorithm. Similarly, the running time per computed coface in Table 4 In Table 5, we present the execution time of the facet generation algorithm applied to simplices of low codimension m, ranging from 1 to 7, in FK-triangulations of high dimensions d (up to d = 400). In Table 6, we present the execution time of the cofacet generation algorithm under the same circumstances. A graphical display of the results of Tables 5 and6 is shown in Figure 10. The total number of intersection points N T = T k ∩ M is thus bounded as follows

N T ≤ K d(d -1)/2 n × i∈I √ dL T ∥u i ∥ D ≤ K n! × d 2 √ dL T R T 2D n . ( 8
)
Here the binomial coefficient arises as the number of choices of n families of hyperplanes.

L C =    √ d+1 2 if d is odd, 1 2
d(d+2) (d+1) if d is even, (9)
and hence L C < √ d+2

2 . We then deduce from ( 8) where x {k} denotes k consecutive coordinates x. This simplex is a simplex in the Coxeter triangulation, as defined in Section A.1.3. In [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] we have seen that the circumcentre of this simplex is c = -d -2i 2(d + 1) , with i ∈ {0, . . . , d}. The circumcentre of a Delaunay simplex is a Voronoi vertex. We recall that All simplices in the star of 0 in the Coxeter triangulation are found by consecutive reflection of a single simplex (in this star) in the hyperplanes of H E C that go through 0, that is the hyperplanes with normals r j,k = e j -e k , with j ̸ = k. See for example [START_REF] Bourbaki | Lie groups and Lie algebras[END_REF][START_REF] Humphreys | Reflection groups and Coxeter groups[END_REF][START_REF] Choudhary | Coxeter triangulations have good quality[END_REF]. We also call these reflections the action of the Weyl group.

N C ≤ K n! × d 2 d(d + 2) 2 √ 2D n N F K ≤ K n! × d 3 √ 2D
The reflection R j,k in a plane that goes through the origin with normal r j,k is given by

R j,k (v) = v -2 v • r j,k r j,k • r j,k r j,k = v -(v • r j,k )r j,k .
We find that

R j,k (c) i = (c -(c • r j,k )r j,k ) i = - d -2i 2(d + 1) - 2j -2k 2(d + 1) (δ ij -δ ik ),

R▶ Definition 1 .S o C G 2 0 2 1 XX: 4 Isomanifold

 114 d and consider the set of hyperplanes H E = {x ∈ R d | ⟨x, u⟩ = k, u ∈ E, k ∈ Z}. Let, in addition, H be the hyperplane of R d+1 of equation ⟨x, 1⟩ = 0 where 1 is the vector of R d+1 whose coordinates are all 1. The Freudenthal-Kuhn triangulation is the hyperplane arrangementH E F K associated to the set of vectors E F K = {e 1 , . . . , e d } ∪ {u i,j = e j -e i | 1 ≤ i < j ≤ d}. The Coxeter triangulation of type Ãd is the hyperplane arrangement H E C in R d+1 associated to the set of vectors E C = {r i,j = e i -e j+1 | 1 ≤ i ≤ j ≤ d} , restricted to H ≃ R d .Two important facts follows. On one hand, because Coxeter and Freudenthal-Kuhn triangulations are related by an affine transformation, they have the same combinatorial structure.We call any triangulation that is the image of a Freudenthal-Kuhn triangulation under an affine transformation a CFK-triangulation. In this paper, we restrict our attention to Coxeter and Freudenthal triangulations since they are the simplest, but any CFK-triangulation could be used. The second fact is that each simplex in such a triangulation can be represented as a cell in an arrangement of d(d -1)/2 families of parallel hyperplanes which are known and do not need to be stored.The next crucial observation relates CFK-triangulations and permutahedra, which allows to represent CFK-triangulations in a compact way (The definition and some combinatorial Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations

Figure 1 ▶ Lemma 3 (▶ Lemma 4 (

 134 Figure 1The permutahedral representation of the simplices in the stars of vertices y and y ′ .

  l denotes the Hessian matrix of second order derivatives, | • | denotes the Euclidean norm of a vector and ∥ • ∥ 2 the operator 2-norm of a matrix.2

Theorem 5 .▶ Theorem 7 .

 57 Assume that the function f has bounded complexity. If the precision of the CFK-triangulation satisfies D = O(d -5/2 ), where the constant in the big O depends on γ max , λ min and α max , then M is a manifold isotopic to the zero set M of f . Moreover, we can bound the Fréchet distance between M and M. The Fréchet distance is a quite strong notion of distance and, in particular, it bounds the Hausdorff distance. ▶ Definition 6 (Fréchet distance for embedded manifolds). Let M a and M b be two homeomorphic, compact submanifolds of R d . Write H for the set of all homeomorphisms from M a to M b . The Fréchet distance between M a and M b is d F (M a , M b ) = inf h∈H sup x∈Ma d(x, h(x)). Assume that the function f has bounded complexity. Then, d F (M, M) = O(D 2 ) where the constant in the big O depends on γ max , λ min and α max .

1 j

 1 generically the intersection of an n-flat H σ with σ. More precisely we note that f | -1 σ (0) is an n-flat if the gradients of f i | σ are linearly independent, which can be easily achieved by perturbing f infinitesimally (or at least its values at the vertices). Let τ d-n j and τ d-n-be faces of σ of dimension d -n and d -n -1. An infinitesimal perturbation of f , can prevent either f | -1 σ (0) from intersecting the faces τ d-n-1 j , or the gradients of f i | σ and the normal spaces of τ d-n j (for each fixed j) from failing to span R d . More precise statements on the geometric and topological stability of the triangulation under perturbations of f can be found in the full version of [13, Section 5]. Because f | -1 σ (0) is (generically) the intersection of an n-flat (H σ )

▶ Proposition 9 .

 9 The time complexity of the algorithm is O (k2 n I|S|) where I is the time complexity of one call of the intersection oracle, and |S| is the number of simplices of dimension k output by the algorithm. Since, the intersection oracle reduces to evaluating f at the k + 1 vertices of σ plus solving a k × k linear system, I = O(k ω ) where ω ≈ 2.375.

▶ Proposition 10 (

 10 Size of the output). Assume that M is contained in the unit cube C d = [0, 1] d , and that any k-flat intersects M at most K times. Writing |S| = N C when T is a Coxeter triangulation and |S| = N F K for a Freudenthal triangulation, we have

▶ Theorem 12 .

 12 where |S| is bounded in Proposition 10. If n = O(1), the combinatorial complexity of M is polynomial in d, and δ = 1/D. We combine Propositions 9, 10, and 11 to obtain our main result. Assume that M is contained in the unit cube [0, 1] d and that any affine k-flat intersects M at most K times (K is usually small, and is in particular a constant for algebraic isomanifolds of bounded degree). Let, in addition, D be the precision required on the approximation (the diameter of a simplex in the ambient triangulation T ). The size of the output, and the time complexity of the algorithm are polynomial in the ambient dimension d, and in δ = 1/D, and exponential in the intrinsic dimension n. The same result holds for the full PL-approximation M of M.

S o C G 2 0 2 1 XX: 10 Isomanifold

 110 approximations are known. Especially important are D-nets[START_REF] Cheng | Manifold reconstruction from point samples[END_REF][START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF]. A D-net consists of a finite number of sample points of M such that no point of M is at distance more than D from a sample point (density condition), and no two sample points are closer than cD for some positive constant c (separation condition). A simple volume argument shows that the size of a D-net of a n-dimensional smooth submanifolds is O(1/D n ) [8, Lemma 5.3]. The sample produced by our algorithm is D-dense on the piecewise linear approximation. This implies that we have a sample that has a Hausdorff distance of D + d F (M, M) to the manifold, where d F (M, M) is bounded in Theorem 7. Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations

3 . 7 Figure 2

 372 Figure 2 The effect of the ambient dimension d and of the intrinsic dimension n on the computation time of the manifold tracing algorithm. The reconstructed manifold in the tests is the n-dimensional sphere embedded in R d . The ambient triangulation used is a Coxeter triangulation of type Ãd . The diameter of the full simplices is fixed for all d.

Figure 3

 3 Figure 3 Comparison of the size of the output of the manifold tracing algorithm using two types of ambient triangulations: a Coxeter triangulation of type Ãd (in blue), and the Freudenthal-Kuhn triangulation of R d (in red) with the same diameter 0.07 √ d of d-dimensional simplices. The reconstructed manifold is the 2-dimensional implicit surface "Chair" embedded in R d given by the equations: (x 2 1 + x 2 2 + x 2 3 -0.8) 2 -0.4 (x3 -1) 2 -2x 2

13 Figure 4

 134 Figure 4 The piecewise-linear approximation of a flat torus embedded in R 10 defined by the equations x 2 1 + x 2 2 = 1, and x 2 3 + x 2 4 = 1, and xi = 0 for i > 4, projected to R 3 . The ambient triangulation used is a Coxeter triangulation of type Ã10 with the diameter of the full-dimensional simplices 0.23. The output size |S| is 509 952. The execution time of the algorithm is 231s. The torus has been rotated and translated in R 10 so that the coordinate axes do not play any special role.

Figure 5

 5 Figure 5 Four views of the flat torus in R 4 given by two equations x 2 1 +x 2 2 = 1, and x 2 3 +x 2 4 = 1 cut by the hypersphere (x1 -1) 2 +x 2 2 +(x3 -1) 2 +x 2 4 = 4, projected to R 3 . The ambient triangulation used is a Coxeter triangulation of type Ã4 with the diameter 0.15 of the full-dimensional simplices. The reconstructed boundary is highlighted in yellow. The size |S| of the piecewise-linear approximation is 14 779. The execution time of the algorithm is 1.84s. The torus has been rotated and translated in R 4 so that the coordinate axes do not play any special role.

  and imaginary part yields a real inhomogeneous polynomial in 2d (real) variables on the domain D i . Taking the real and imaginary parts of the function yields two real functions f R,i and f I,i on D i . As real sets, the projective variety f = 0 on CP d and the intersection of the sets f R,i = 0 and f I,i = 0 on D i for each i (with the boundaries identified) are the same. We can therefore apply the tracing algorithm to each isomanifold (f R,i = 0, f I,i = 0) of D i independently. Since their boundaries coincide, we can then glue these isomanifolds along their boundary to obtain a PL-approximation of the projective variety f = 0. This, for example, allows to recover the Euler characteristic of f = 0 on CP d .

▶ Definition 15 (S o C G 2 0 2 1 XX: 22 Isomanifold

 15122 Ordered partition). Let T be a finite non-empty set, |T | its cardinality, and l ≤ |T | a positive integer. An ordered partition of T in l parts is a collection of l indexed subsets ω = (ω 1 , . . . , ω l ), such that ω i ⊆ T , and {ω 1 , . . . , ω l } is a partition of T . The ω i are called the parts and are ordered by their index. We write OP l [d] for the set of ordered partitions of [d] with l parts, and just OP [d] for the set of all ordered partitions of [d]. Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations

  where |ω p | denotes the length of the p-th part of the ordered partition, and P(n) the permutahedron of dimension n -1. Proof of Corollary 19. Write l = d -i. Since the number of vertices of the product of two polytopes is the product of the vertices, and a (n -1)-dimensional permutadedron has n! vertices, we see that the total number of vertices of a i-face of a d-dimensional permutahedron corresponding to an ordered partition ω = (ω 1 , . . . , ω l+1 ) is l+1 p=1 (|ω p |!). Let 1 ≤ j < k ≤ d be integers such that j + k = d + 1. By definition j!k! < (j -1)!(k + 1)!, and thus j!k! ≤ 1!d!. Generalizing this, we see that the product of the |ω p |! is maximal when all parts are singletons except the biggest part which has d + 1 -l elements. Therefore l+1 p=1 (|ω p |!) ≤ (d -l + 1)!.

and thus d+1 2 ≥(▶

 2 l + 1, then we see thatl+1 p=1 (|ω p |!) ≥ 2 l+1 . If ⌊ d+1 l+1 ⌋ = 1, we have that d+1 2 < l + 1 ≤ d + 1, and thus r ′ = d + 1 mod l + 1 = d -l. |ω p |!) ≥ 2 d-l .Because l + 1 = d -l, or 2l + 1 = d is precisely the point where you go from the first to the second case we see that l+1 p=1 (|ω p |!) ≥ 2 min{l+1,d-l} ◀ Corollary 21. The number of facets of an i-face σ of a d-permutahedron is at most 2 i+1 .

  the previous two corollaries. For i = 0, the bound in Corollary 22 is the same as the upper bound in Corollary 19. For i = j -1, the bound is comparable but weaker than the bound in Corollary 21.Proof of Corollary 22.Let σ be a j-face of the d-permutahedron. Write F i,σ for the set of i-faces of σ, and c v for the number of i-cofaces of a vertex v of σ, i.e. the number of simplices of F i,σ that contain v. For τ ∈ F i,σ , we write p τ for the number of vertices of τ . By double counting the incidences between vertices, and i-faces inside σ, we haveτ ∈Fi,σ p τ = v∈F0,σ c v .Now observe that the d-permutahedron is a simple polytope (this follows from the fact that its dual is simplicial since it is a star in the FK-triangulation). The faces of simple polytopes are also simple polytopes which implies that the vertices of a j-face are incident to j i faces of dimension i [14, Lemma 7.1.14]. Moreover, |F 0,σ | ≤ (j + 1)! by Corollary 19. Hence v∈F0,σ

Figure 8

 8 Figure 8The ordered partitions associated to the faces of the FK-triangulation of R 2 that have the same minimal vertex v0 (circled).

S o C G 2 0 2 1 XX: 28 Isomanifold

 128 Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations We will call any triangulation of R d that is the image of a Freudenthal-Kuhn triangulation under a non-degenerate affine map a Coxeter-Freudenthal-Kuhn triangulation, or CFKtriangulation for short. This includes the Coxeter triangulation of type Ãd (as embedded in R d ).

Figure 9

 9 Figure 9The Coxeter and Freudenthal-Kuhn triangulations in the plane.

[ 14 ,▶ Proposition 33 .S o C G 2 0 2 1 XX: 30 Isomanifold▶ Proposition 34 .

 143313034 Section 11.3]). The star of a vertex in a CFK-triangulation is combinatorially dual to a permutahedron.This proposition explains the nomenclature permutahedral representation. This is equivalent to the following more geometric result: Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations The Voronoi cell of a Coxeter triangulation of type Ãd is a permutahedron.

Point location.

  Given a point x ∈ R d Lemma 24 tells us how to find the canonical permutahedral representation of the simplex in which x is contained. The complexity of point location is dominated by the sorting of the z i = x i -⌊x i ⌋, which takes O(d log d) time, and requires O(d) space.

▶ Lemma 35 (

 35 takes time O(k + s) where s = l+1 k+1 is the number of subsets. Computing, for each such k-simplex its permutahedral representation takes O(d) time. Face computation (full version)). Let σ be an l-simplex in the FK-triangulation of R d given by its canonical permutahedral representation. Computing the canonical permutahedral representations of all its k-faces can be done in time O(ds), where s = l+1 k+1 is the number of k-faces of an l simplex. The space complexity of the algorithm is O(l).

1 . 2 .

 12 The geometry of the triangulation is specified by the affine transformation that maps the FK-triangulation of R d to the CFK-triangulation. The affine transformation is given by a d × d matrix Λ, and a d-vector b. For the FK-triangulation, Λ is the identity matrix, and b = 0; therefore no storage is required. For the Coxeter triangulation of type Ãd , Λ is sparse, as can be seen by inspection of the proof of Lemma 28.

S o C G 2 0 2 1 XX: 34 Isomanifold

 134 Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations The next results are motivated by the problem of tracing a manifold of low dimension m embedded in R d for high d. The crucial operations in this context consist in computing the facets and cofacets of simplices of codimension m in a triangulation of R d , as is clear from Algorithm 1 .

Figure 10 9 .

 109 Figure 10Graphical display of the results of Tables5 and 6.

2 .

 2 Consider now more specifically Coxeter triangulations of type Ãd and F K-triangulations. It follows from Section A.1 that R C = R F K = √ The longest edge L F K in a Freudenthal-Kuhn triangulation has length at most (in fact exactly) √ d since each simplex is contained in a cubical cell of the d-dimensional unit grid. Furthermore, it is proved in [18, point 6 of Ãd in Section 6] that the longest edge length in the Coxeter triangulation of type Ãd is

1 XX: 36 IsomanifoldC Alternative Proof of Proposition 34 Proof of Proposition 34 .

 1363434 Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations Proof of Proposition 11. Let σ be a k-simplex of a CFK-triangulation that intersects M, and let σ * be its dual cell. By definition, σ * is a m-dimensional face of T * , the polytopal cell complex dual to T . The collection of all σ * associated to the k-simplices σ of T that intersect M form a cell complex M * dual to M. To bound the number of faces of all dimensions of the PL-approximation M, it is therefore sufficient to bound the number of faces of M * .Each d-dimensional cell in T * is a permutahedron (Proposition 33). Hence, σ * is a n-face of a d-permutahedron. The number of faces of σ * of dimensions 0 to n -1 (or equivalently the number of cofaces of σ of dimensions n + 1 to d) is: p i,j denotes the number of i-faces of a j-face of the d-permutahedron and is bounded in Corollary 22. The last equality can be easily verified using Mathematica. The overall combinatorial complexity of M is therefore|S| × 3 n -1 2 n (n + 1)!,where S is bounded in Proposition 10. ◀ We start by recalling a number of results. Let P = {(x i ) ∈ R d+1 | i x i = 0} and consider the d-simplex with vertices u k in P .

G 2 0 2 1 XX:8 Isomanifold Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations seed

  Given a k-simplex σ of T , decide whether σ intersects H σ and, in the affirmative, report the corresponding vertex σ ∩ H σ = σ ∩ M. simplices. If M consists of multiple connected components, then a seed simplex must be provided per each connected component and we proceed in the same manner for each component. So we will assume for now that M is connected.

		Algorithm 1 Manifold tracing algorithm
	218 input : the permutahedral representation of a triangulation T of R d ,
		219 the codimension of the isomanifold k = d -n,
		220 a seed k-simplex τ 0 that intersects	M
	221 oracle : 222 output : Set S of the simplices in T of dimension k that intersect M, represented by
		their permutahedral representation, and the corresponding set M0 of
		intersection points
	223 Initialize the queue Q and the set S with τ 0
	224 while the queue Q is not empty do
	225	Pop a k-dimensional simplex τ from Q
	226	foreach cofacet ϕ of τ do
	227	foreach facet σ of ϕ do
	228	if σ does not lie in S and intersects M (which can be decided using the
		oracle) then
	229	Insert σ into the queue Q
	230	Insert σ into S together with the intersection point provided by the
		oracle
	The algorithm essentially computes the set S of k-simplices of T that intersect M. The
	elements of S are in 1-1 correspondence with the vertices of M thanks to the Genericity
	hypothesis. The so-called intersection oracle is a basic ingredient of the algorithm:
	Intersection oracle: Given a k-simplex σ of T , decide whether σ intersects H σ and, in the
	affirmative, report the corresponding vertex σ ∩ H σ .
	It is easy to see that the intersection oracle reduces to solving a linear system. Indeed,
	generically, a vertex is the intersection of a k-simplex σ of T with the m-flat H σ that
	interpolates f inside σ. One can compute the barycentric coordinates of σ ∩ H σ by solving
	a linear system of k equations, and k unknowns. It then remains to check whether the
	barycentric coordinates are all non-negative (to ensure that the intersection point lies inside
	σ). It follows that the intersection oracle reduces to evaluating f at the k + 1 vertices of σ
	plus solving a k × k linear system.
	In addition, we need to provide a set of k-simplices of T to initialize the tracing. These
	simplices must intersect all the connected components of the isomanifold and are called
			S o C

  k+1 of (k + 1)-faces of the triangulation T that intersect ∂ M. The computation of S k+1 is done by the following simple modification of Algorithm 1: if the k-dimensional facet σ of τ intersects f -1 (0) at a point x such that f∂ | τ (x) < 0 (i.e. x is not in M), we then compute the intersection point of τ with f -1 ∂ (0), and put τ in S k+1 .As for the case of manifolds without boundary (see the discussion at the end of Section 3.2), the algorithm traverses (and therefore computes) the 1-skeleton of M. Under the Genericity Hypothesis 8, the vertices of M1 are in bijection with the simplices of S k ∪ S k+1 . The edges are obtained by applying the following rules below (we identify a simplex in S k (resp. S k+1 )and the intersection point S k ∩ M (resp. S k+1 ∩ ∂ M):1. Two simplices σ 1 ,and σ 2 of S k are joined by an edge in M1 if and only if there exists a simplex in T k+1 with faces σ 1 and σ 2 . 2. Two simplices τ 1 , and τ 2 of S k+1 are joined by an edge in ∂M 1 if and only if there exists a simplex in T k+2 with faces τ 1 and τ 2 . 3. A simplex σ of S k , and a simplex τ of S k+1 are joined by an edge in ∂M 1 if and only if

J.-D. Boissonnat, S. Kachanovich, and M. Wintraecken XX:11 S σ is a facet of τ .

XX:24 Isomanifold Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations Proof.

  Write l = d -i. We first recall a set of d > 2 objects can be subdivided in two non-empty ordered subsets A and B in 2 d -2 ways. This is not hard to see. Because we pick for each element if it will be put in A or B there are 2 d possibilities. Excluding that A or B is empty gives 2 d -2. Let ω = (ω 1 , . . . , ω l ) again be an ordered partition. To find a refinement of ω in l + 1 parts, we need to first pick a 1 ≤ p ≤ l, such that |ω p | > 1, and then we need to break ω p up into two (ordered) parts, for which there are 2 |ωp| -2 possibilities as we have seen above. This means that ifI = {p | 1 ≤ p ≤ l, |ω p | > 1}, the number of ≤ s < t ≤ d be integers such that s + t = d + 1. Then 2 s + 2 t < 2 s-1 + 2 t+1 .Generalizing this, we see that the sum of the 2 |ωp| -2 is maximal when all |ω p | = 1 except the biggest part which has d -l + 1 = i + 1 elements.

	refinements is
	2 |ωp| -2.
	p∈I
	Let now 1

◀ ▶ Corollary 22. Let p i,j denote the number of i-faces of a j-face of the d-permutahedron.

  to obtain the canonical permutahedral representation of the coface. The next Lemma follows. The bound on the number s of cofaces follows, by duality, from Corollary 22. ▶ Lemma 36 (Coface computation (full version)). Let σ be a k-simplex in the FK-triangulation of R d given by its permutahedral representation. Computing the permutahedral representations of all its l-cofaces can be done in time O(ds), where

Table 2

 2 shows the same running times per computed face. As we can see, except for the case l = k, the running time per computed face is around 2µs.

S o C G 2

0 2 1 XX:32 Isomanifold Tracing in R d , using Coxeter-Freudenthal-Kuhn Triangulations

	Face dimension k 22	23	24	25	26	27	28	29
		22	0.006
	Simplex dimension l	23 24 25 26 27 28 29	0.042 0.503 4.88 33.76 162.114 35.108 0.006 0.05 0.645 5.697 885.293 190.441 40.856 0.008 0.058 0.697 6.824 3420.99 973.455 246.88	0.008 0.062 0.758 6.906 49.896	0.008 0.064 0.739 6.657	0.008 0.058 0.006 0.735 0.058 0.006
		30	11904.7 4175.92 1247.97 275.776 50.932 7.348 0.778 0.058

Table 1

 1 Total running time of the face generation algorithm (in milliseconds).

	Face dimension k 22	23	24	25	26	27	28	29
		22	0.006
	Simplex dimension l	23 24 25 26 27 28 29	0.0018 0.006 0.0017 0.002 0.0019 0.002 0.0019 0.0019 0.002 0.008 0.0022 0.008 0.0023 0.008 0.0016 0.0017 0.0021 0.002 0.0023 0.006 0.0019 0.0016 0.0017 0.0019 0.0018 0.002 0.0017 0.0016 0.0017 0.0018 0.0016 0.0017 0.0019 0.006 0.006
		30	0.0015 0.0016 0.0017 0.0016 0.0016 0.0016 0.0017 0.0019

Table 2

 2 Running time of the face generation algorithm per computed face (in milliseconds).

Table 3

 3 is around 2µs with the exception of when k is close to l. Total running time of the coface generation algorithm (in milliseconds).

	Coface dimension l 23	24	25	26	27	28	29	30
		22	0.11	1.274 9.577 43.848 86.699 96.407 59.935 15.487
	Simplex dimension k	23 24 25 26 27 28 29	0.043 0.114 0.729 3.499 0.047 0.1 0.381 0.046 0.097 0.047	9.337 1.183 0.23 0.076 0.049	13.523 10.058 3.049 2.132 1.871 0.653 0.423 0.426 0.193 0.128 0.15 0.093 0.069 0.081 0.063 0.047 0.061 0.054 0.05 0.053
		30							0.05
	Coface dimension l 23	24	25	26	27	28	29	30
		22	0.002 0.0013 0.0013 0.0015 0.0016 0.0016 0.0016 0.0017
	Simplex dimension k	23 24 25 26 27 28 29	0.042 0.003 0.045	0.0017 0.0016 0.0016 0.0016 0.0016 0.0017 0.004 0.0019 0.0017 0.0017 0.0017 0.0018 0.045 0.0053 0.0025 0.002 0.0019 0.0022 0.047 0.0073 0.0035 0.0028 0.0036 0.048 0.0103 0.0058 0.0068 0.048 0.0145 0.0133 0.05 0.026
		30							0.05

Table 4

 4 Running time of the coface generation algorithm per computed face (in milliseconds).

Table 5

 5 Average running times in milliseconds of the facet generation algorithm.

	Ambient dimension d 50	100	150	200	250	300	350	400
	Face codimension	m	1 2 3 4 5 6 7	0.166 0.612 1.438 2.862 5.376 8.69 0.166 0.643 1.417 2.858 5.607 8.375 11.806 16.261 12.184 15.924 0.168 0.607 1.395 2.888 5.866 8.232 12.008 16.527 0.162 0.589 1.373 2.864 5.491 8.447 11.936 16.08 0.154 0.587 1.349 2.76 5.77 8.371 11.814 15.88 0.148 0.579 1.321 2.737 5.735 8.351 12.038 15.798 0.136 0.575 1.313 2.553 5.701 8.313 12.11 15.754
	Ambient dimension d 50	100	150	200	250	300	350	400
			1	0.068 0.134 0.228 0.281 0.423 0.605 0.611 0.848
	Simplex	codimension m	2 3 4 5 6	0.082 0.17 0.098 0.194 0.303 0.401 0.525 0.733 0.866 1.124 0.267 0.341 0.483 0.723 0.731 0.966 0.112 0.226 0.351 0.467 0.665 0.806 0.974 1.295 0.132 0.265 0.423 0.545 0.966 0.928 1.128 1.477 0.162 0.329 0.515 0.713 0.948 1.124 1.361 1.76
			7	0.2	0.415 0.651 0.878 1.166 1.421 1.784 2.283

Table 6

 6 Average running times in milliseconds of the cofacet generation algorithm.

Acknowledgements

We thank Dominique Attali, Guilherme de Fonseca, Arijit Ghosh, Vincent Pilaud and Aurélien Alvarez for their comments and suggestions. We also acknowledge the reviewers.

Funding The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. 339025 GUDHI (Algorithmic Foundations of Geometry Understanding in Higher Dimensions). Jean-Daniel Boissonnat: Supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. Mathijs Wintraecken: Supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 754411.

which permutes the jth and kth coordinate of c. Here we used the upper index i to denote the ith coordinate. Using the cycle notation for the permutation group, see for example [START_REF] Armstrong | Groups and Symmetry[END_REF]Chapter 6], this coincides the 2-cycle (j k). Let now

with {π i } some permutation of {0, . . . , d}. We find that

which again permutes the jth and kth coordinate. Now recall that all permutations are generated by 2-cycles, see for example [START_REF] Armstrong | Groups and Symmetry[END_REF]Theorem 6.1]. This implies that, for any permutation π, we can find c π from c by the action of the Weyl group. This also means that we have explicitly described the Voronoi cell of 0 in the Coxeter triangulation of type Ãd as a permutahedron. Because of symmetry, this now holds for any Voronoi cell. ◀