Jean-Daniel Boissonnat
email: jean-daniel.boissonnat@inria.fr

Siargey Kachanovich
email: siargey.kachanovich@inria.fr

Mathijs Wintraecken
email: m.h.m.j.wintraecken@gmail.com

{1

Sampling and Meshing Submanifolds in High Dimension

Keywords: Subject Classification Theory of computation → Computational geometry Coxeter triangulation, Kuhn triangulation, permutahedron, PL-approximations, isomanifolds { 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}

This paper presents a rather simple tracing algorithm to sample and mesh an m-dimensional submanifold of R d for arbitrary m and d. We extend the work of Dobkin et al. to submanifolds of arbitrary dimension and codimension. The algorithm is practical and has been thoroughly investigated from both theoretical and experimental perspectives. The paper provides a full description and analysis of the data structure and of the tracing algorithm. The main contributions are :

1. We unify and complement the knowledge about Coxeter and Freudenthal-Kuhn triangulations.

2. We introduce an elegant and compact data structure to store Coxeter or Freudenthal-Kuhn triangulations and describe output sensitive algorithms to compute faces and cofaces or any simplex in the triangulation. 3. We present a manifold tracing algorithm based on the above data structure. We provide a detailled complexity analysis along with experimental results that show that the algorithm can handle cases that are far ahead of the state-of-the-art.

Introduction

This paper presents a rather simple algorithm to sample and mesh an m-dimensional submanifold of R d for arbitrary m and d. This fundamental problem finds applications in various fields like numerical analysis (to solve nonlinear differential equations) [START_REF] Ambrosio | Level set approach to mean curvature flow in arbitrary codimension[END_REF], dynamical systems (to approximate invariant manifolds) [START_REF] Michael | The computation of fixed points and applications[END_REF][START_REF] Haller | Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction[END_REF]), chemistry (to study the energy landscape of molecules) [START_REF] Martin | Topology of cyclo-octane energy landscape[END_REF], robotics (to describe the configuration space of mechanical systems) [START_REF] Lavalle | Planning Algorithms[END_REF], computer vision and graphics (to visualize time-varying and higher dimensional data) [START_REF] Bhaniramka | Isosurfacing in higher dimensions[END_REF][START_REF] Min | Simplicial isosurfacing in arbitrary dimension and codimension[END_REF].

State-of-the-art. The problem of triangulating differentiable manifolds has a long history in Mathematics dating back to the work of Cairns [START_REF] Cairns | On the triangulation of regular loci[END_REF], Whitehead [START_REF] Whitehead | On C 1 -complexes[END_REF] and Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF].

More recently, the problem has received a lot of attention for surfaces of R 3 in the Computational Geometry and Computer Graphics communities. Among the widely used methods are the Marching cube algorithm [42] and Delaunay refinement [START_REF] Cheng | Delaunay Mesh Generation. Chapman and Hall / CRC computer and information science series[END_REF]. In higher dimensions, some early work 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]. 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 Coxeter triangulations and their potential use for contour tracing. Although the authors only considered the case of curves (m = 1), it was claimed that the method could be "immediately generalized" to submanifolds of higher codimensions. However there has been only few works in that direction. This situation might be explained by the fact that extending the Marching cube algorithm to higher dimensions seems infeasible due to the large number of configurations that should be stored in a lookup table, and that no efficient data structure was known to store and query triangulations in high dimensions. The most recent work we are aware of is the work of Bhaniramka et al. [START_REF] Bhaniramka | Isosurfacing in higher dimensions[END_REF] which is limited to hypersurfaces (m = d -1) and the work of Min [START_REF] Min | Simplicial isosurfacing in arbitrary dimension and codimension[END_REF]. Min's method applies to submanifolds of any dimension and codimension. It uses an ambient triangulation instead of a cubical grid (the same Freudenthal-Kuhn triangulation used in this paper). However, in their analysis, they consider d as a constant and only report experimental results in dimensions 3 and 4.

The problem we consider is also related to the problem of manifold reconstruction from point samples [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF][START_REF] Boissonnat | Manifold reconstruction in arbitrary dimensions using witness complexes[END_REF][START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF]. A major difference however is that in manifold reconstruction, a sample is given as input while here we have to construct the sample using an oracle that queries the manifold. Manifold sampling is another related problem which is of fundamental algorithmic significance in statistics. Yet, not much is known beyond the convex case and the case of hypersurfaces (m = d -1) [20,[START_REF] Narayanan | Sampling hypersurfaces through diffusion[END_REF].

Contributions. This paper is the first of a series of related papers to fill the gap. In this paper, we extend the work of Dobkin et al. [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF] to submanifolds of arbitrary codimension.

The algorithm is practical and has been thoroughly investigated from theoretical and experimental perspectives. This paper provides a full description and analysis of the data structure and of the tracing algorithm.

Guarantees on output of the the algorithm are established in two companion papers, one for the case of isomanifolds [START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF] and one for general smooth submanifolds [START_REF] Boissonnat | Triangulating submanifolds: An elementary and quantified version of Whitney's method[END_REF]. Specifically, for isomanifolds we prove [START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF] that the output is a PL-manifold that has the same topology type as M, and whose Fréchet distance to M is small. Implementation details and experimental results will be discussed in [START_REF] Boissonnat | Tracing isomanifolds of arbitrary dimension and codimension[END_REF]. The case of submanifolds with boundaries and, more generally, of stratified manifolds can be handled in very much the same way [START_REF] Boissonnat | Tracing isomanifolds of arbitrary dimension and codimension[END_REF][START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF].

The content of this paper is as follows. In Section 2, we first discuss Freudenthal-Kuhn and Coxeter triangulations, the latter exclusively of type Ãd . triangulations are the same up to a linear transformation, as noted by [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF] and fully proved in the appendix. This allows us to use on one hand the nice geometric properties of Coxeter triangulations of type Ãd , where each simplex is very well shaped (large volume compared to longest edge length) and all simplices are identical up to reflections, and on the other hand the simple combinatorial definitions of the Freudenthal-Kuhn triangulation. Although most ideas in this section were known prior to this work, they were disseminated in many different areas and difficult to access. We elucidate those ideas, provide full proofs and combine them so as to extend them to arbitrary dimensions when necessary.

We then introduce a new data structure to compactly store Coxeter or Freudenthal-Kuhn triangulations. The data structure is an elegant and efficient representation of the combinatorial structure of those triangulations. We present a point location algorithm and describe output sensitive algorithms for computing faces and cofaces of simplices of all dimensions in the triangulation.

In Section 3, we present a submanifold tracing algorithm based on the above data structure. The algorithm works for smooth submanifolds of any dimension and codimension.

Starting from a given seed point, the algorithm probes the manifold using an oracle and outputs both a sample and a PL approximation of the manifold. A distinctive property of our algorithm, when compared to previous methods that work for submanifolds of arbitrary dimension and codimension [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF][START_REF] Min | Simplicial isosurfacing in arbitrary dimension and codimension[END_REF], is that its complexity depends mostly on the intrinsic dimension of the manifold (see Theorem 24 for a precise statement). Furthermore, using dimensionality reduction techniques, we can completely remove the dependency of the result on the ambient dimension (Theorem 26). The algorithm is quite simple and our implementation can handle cases that are far beyond what was possible before (Section 4).

2

Coxeter-Freudhental-Kuhn triangulations

Permutahedra

We write [k] = {1, . . . , k} and [k, l] = {k, . . . , l}.

Definition 1 (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].

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]. For example, let T = {1, 2, 3, 4, 5}. If {1, 5}, {2, 3, 4} is an ordered partition of T , then {1, 5}, {2}, {4}, {3} is a refinement of {1, 5}, {2, 3, 4}. Definition 1.3.10. We denote by OP (T ) the set of all ordered partitions of a finite set T .

The refinement partial order v defines a partial order on the set of ordered partitions OP (T ) [MK92, Section 3.1.2]. Consequentially, (OP (T ), v) is a partially-ordered set, called refinement poset.

Now, we are ready to state the main result of this section.

33 124

We introduce the convention that z d+1 = 0. We associate to x the ordered partition ω =

125

(ω 1 , . . . , ω l+1 ) of [d + 1] where the ω i are obtained by sorting the z i in decreasing order.

126

Specifically, with

ω i = {ω i (1), . . . , ω i (m i )}, we have 127 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. ( 1
)
Lemma 9. Suppose that ω = (ω 1 , . . . , ω l+1 ) is an ordered partition of [d + 1] and let σ = {v 0 , . . . , v l } be the l-simplex whose vertices are the points

v 0 = ( x 1 , . . . , x d ), v i = v i-1 + i∈ωi e i i = 1, . . . , l. ( 2
)
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 and d + 1 ∈ ω l+1 ), satisfy [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF].

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 , (3)
for some λ i > 0 satisfying l i=0 λ i = 1. We write

α ω l (1) = • • • = α ω l (m l ) = λ l . . . α ω1(1) = • • • = α ω1(m1) = λ l + • • • + λ 1 (4)
By construction α ω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 (3), for all i. Moreover, because

λ l + • • • + λ i > λ l + • • • + λ i-1 , we see that (1) is satisfied.
Conversely, given a point z such that its coordinates satisfy (1), we can read of its barycentric coordinates with respect to the v i , as defined by (2), from (4).

Theorem 10. 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.

Proof. Lemma 9 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 σ. Because x has uniquely defined barycentric coordinates with respect to the vertices of σ it is mapped to a unique point in σ.

Hence the partition we have defined is a well-defined triangulation of R d .

Remark. We note that, by construction, v 0 in Lemma 9 is the smallest vertex of σ in the lexicographical order. Lemma 9 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.

S o C G 2 0 2 0 00:6

Sampling and Meshing Submanifolds

Cycles and the permutahedron. This monotone path can be made into a cycle using the extra vector e d+1 , introduced by Eaves [START_REF] Eaves | A course in triangulations for solving equations with deformations[END_REF], because by construction

v 0 = v l + i∈ω l+1 e i ,
with ω as in Definition 8. 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 11 (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. ( 5
)
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. In Lemma 14 and Proposition 15 we'll 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 12 (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 . ( 6
)
Here we use the convention that the sum from 1 to 0 is empty. We write (v 0 , ω ) = (v 0 , ω)⊕k.

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

Proof. Follows by inserting (6) in [START_REF] Bhaniramka | Isosurfacing in higher 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 14. 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 12 and Lemma 13, (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 9 and Theorem 10.

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 (5) 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 . For example, when looking at the two dimensional example in Figure 2, we see that the edges that contain y in the red triangle with Because the combinatorial structure of the faces is compatible with the permutahedron, Lemma 14 immediately gives:

Proposition 15. The star of v 0 is dual to a permutahedron (combinatorially).

This proposition also explains the nomenclature permutahedral representation. τ is a k-face of σ. We then compute the canonical permutahedral representation of all those k-faces τ .

Basic operations

Point location

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 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 14 to obtain the canonical permutahedral representation of the coface.

Lemma 17.

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 s is the number of l-cofaces of a k-simplex in the FK-triangulation.

The space complexity of the algorithm is O(d).

Coxeter triangulations of type Ãd

The Freudenthal-Kuhn triangulation is closely related to the Coxeter triangulation [START_REF] Harold | Discrete groups generated by reflections[END_REF] of type Ãd . There are many equivalent ways to define the Coxeter triangulation of type Ãd , see [START_REF] Bourbaki | Lie groups and Lie algebras[END_REF][START_REF] Choudhary | Coxeter triangulations have good quality[END_REF][START_REF] Humphreys | Reflection groups and Coxeter groups[END_REF]. We recall the following:

u 0 = 0 {d+1} u k = - d + 1 -k d + 1 {k} , k d + 1 {d+1-k} , k ∈ [d],
where x {k} denotes k consecutive coordinates x. The Coxeter triangulation of type Ãd in P is found by consecutively reflecting the simplex in its faces.

The following lemma relates Coxeter and FK-triangulations and was first stated in [START_REF] Dobkin | Contour tracing by piecewise linear approximations[END_REF].

Lemma 19. The Freudenthal-Kuhn triangulation and Coxeter triangulation of type Ãd

are identical up to a linear transformation.

A proof can be found in Appendix A. We will now call any triangulation of Euclidean space that is the image of a Coxeter triangulation under a non-degenerate affine map a Coxeter-Freudenthal-Kuhn triangulation, or CFK-triangulation for short. Moreover, we have Lemma 20 ([15]). The Coxeter triangulation of type Ãd is a Delaunay triangulation.

We note that Proposition 15, together with Lemmas 19 and 20, give an alternative selfcontained proof of the known fact [17, Chapter 21, Section 3.F] (also proved in Appendix A for completeness) that the Voronoi cell of a vertex in Ãd is a (combinatorial) permutahedron.

The simplices in the Coxeter triangulation of type Ãd have extremely good quality [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF].

For example, the volume compared to the longest edge length to the d-th power is large. As we will see, the exceptional quality of Coxeter improves the running of our algorithms.

Data structure for storing CFK-triangulations

To store an ambient CFK-triangulation for the manifold tracing algorithm in Section 3, we use the following data structure. This data structure contains information on both the combinatorial structure and the geometry of the triangulation. The combinatorics of the triangulation is given through the canonical permutahedral representation of its simplices and the algorithms from Section 2.2. 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.

Remark. Matrix Λ is used to compute the coordinates of the vertices of simplices and, for the most useful cases in practice, needs not to be explicitly stored. For the FK-triangulation, Λ is the identity matrix and b = 0, therefore no storage is required. For the Coxeter triangulation of type Ãd , we can directly access the coordinates of vertices as given in Definition 18.

Sampling and meshing submanifolds

In this section, we describe an algorithm that will compute a PL-approximation of an msubmanifold of R d for arbitrary d and m ≤ d. The algorithm can be considered as an alternative to the Marching Cube algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] where the usual cubical grid is replaced by a CFK -preferably the Coxeter -triangulation of the ambient space. Taking a 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]. By taking as a triangulation of the ambient space a CFK-triangulation, we keep two main advantages of using grids: very limited storage and fast basic operations. Insert ρ together with the intersection point to the output set S

Manifold tracing algorithm

Let M be an m-dimensional compact submanifold of the Euclidean space R d . Both m and d are known but arbitrary and will be considered as parameters in the complexity analysis. The algorithm will use a CFK-triangulation T of R d , which is stored using the data structure from Section 2.5. We assume that the manifold M and the triangulation T satisfy a genericity hypothesis:

Hypothesis 21 (Genericity). The manifold has an empty intersection with all simplices of dimensions strictly lower than k in the triangulation T . The intersection of the manifold

M and any k-dimensional simplex in the triangulation T is a single point.

Remark. It turns out [START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF] that for an isomanifold f -1 (0) it suffices to find the intersection points of f -1 PL (0) with the k-simplices under very weak conditions. Here f PL denotes the function that is linear on every simplex in T and coincides with f on the vertices of T . We stress that f -1 PL (0) satisfies the genericity hypothesis with probability one.

We assume that we know a point on the manifold x 0 ∈ M, from which the algorithm starts. If M consists of multiple connected components, then a seed point per each connected component must be provided and we proceed in the same manner for each component. So we will assume for now that M is connected.

In addition, we assume that the manifold M can be accessed through an oracle that

allows us to answer whether a k-simplex in the triangulation T intersects the manifold M.

Here, k = d -m is the codimension of M. In the following, we will refer to this oracle as the intersection oracle.

The algorithm is described as Algorithm 1. We first translate the coordinate frame so that x 0 is the barycenter of a k-simplex of T (any such simplex is fine). This simplex is put in the set S of the simplices in T of dimension k = d-m that intersect M. Then, given such a simplex, we look at all its cofacets that have not been considered yet and consider all the facets of those cofacets that have not been considered yet. This can be done using a queue Q of simplices to consider. Each of these simplices is queried with the intersection oracle and, if it is found to intersect M, it is added to S. Upon termination, S contains all the k-dimensional simplices of T that intersect M. Since, by our genericity assumption, each k-simplex in S intersects M in a single point, |S| is also the size of the sample produced by our algorithm. A better approximation of the sample is of course possible if we have at our A polyhedron can be deduced from S by taking the dual faces of the simplices in S. A more precise approximation can be obtained if, in addition to the intersection oracle, we can also compute the intersection points S ∩ M. This will be described in full detail in a companion paper.

Complexity analysis

We can easily bound the complexity of the manifold tracing algorithm as a function of the size of the output.

O(d2 m ) + O(dk2 m ) + O(k2 m I) = O(k2 m (d + I)) = O(k2 m I).
Since there are |S| iterations of the while loop, the result follows.

We will now express the size of the output in terms of quantities that depend on the manifold and the resolution of the triangulation.

Proposition 23 (Size of the output).

|S| ≤ C Θ √ m 2πe k k/2 2 δ m vol m (M) = O 2 O(d log d) vol m (M) δ m .
where:

vol m (M) is the m-dimensional volume of M,
δ is the diameter of the d-simplices of T and a measure of the resolution of T ,

V is the volume of any d-simplex of T , and

Θ = V δ d its fatness,
C is a constant that does not depend on d, m or δ.

Proof. Let N be the set of the d-dimensional cofaces of the simplices in S, and let N be the cardinality of N . In the proof we will use constants C 1 , C 2 , C 3 that are constants that do not depend on d, m nor δ.

Upper bound on N . Write M δ for the tubular neighbourhood of M of radius δ, i.e. the set of points at distance at most δ from M. Since the d-dimensional simplices in N have pairwise disjoint interiors and all lie inside M δ , we have

N • V ≤ vol d (M δ ). (7)
S o C G 2 0 2 0 00:12

Sampling and Meshing Submanifolds

According to the tube formula of Weyl [START_REF] Weyl | On the volume of tubes[END_REF][START_REF] Gray | Comparison theorems for the volumes of tubes as generalizations of the weyl tube formula[END_REF] and writing B k for the volume of the unit ball of dimension k, there exists constants C 1 and C 2 such that

vol d (M δ ) ≤ C 1 B k δ k vol m (M) ≤ C 2 2πe k k/2 δ k vol m (M) (8)
By combining the two inequalities ( 7) and ( 8), we get:

N ≤ vol d (M δ ) V ≤ C 2 2πe k k/2 δ k V vol m (M) ≤ C 2 2πe k k/2 1 Θ vol m (M) δ m . ( 9
)
where Θ is the fatness of the simplices in the triangulation T . Note that the dependency of

N
The number of d-cofaces of a k-simplex is given by Corollary 6 applied to the dual Voronoi face of τ . Hence there exists a constant C 3 such that for any k-dimensional simplex τ ∈ S, we have:

|cof(τ, d)| ≥ C 3 d m + 1 m+1 .
On the other hand, for each d-dimensional simplex σ ∈ N , we have :

|fac(σ, k) ∩ S| ≤ |fac(σ, k)| = d + 1 m + 1 .
Equation [START_REF] Boissonnat | The topological correctness of PLapproximations of isomanifolds[END_REF] becomes

C 3 d m + 1 m+1 |S| ≤ τ ∈S |cof(τ, d)| = σ∈N |fac(σ, d) ∩ S| ≤ d + 1 m + 1 N,
from which we get using ( 9) that there exists a constant C such that

|S| ≤ C 2πe k k/2 1 Θ vol m (M) m -1/2 (2/δ) m .
Bound on fatness. The fatness term Θ in the expression in Proposition 23 depends on the choice of the triangulation T . The fatness Θ CT of the d-dimensional simplices in the Coxeter triangulation is given by (see [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF])

1 Θ CT = O d (d+1)/2 d! 2 d = O(2 O(d log d) ). ( 11
)
while the fatness Θ F KT of the d-dimensional simplices in the FK-triangulation of R d is given by (see [START_REF] Dahmen | On the linear independence of multivariate B-splines, I. Triangulations of simploids[END_REF])

1 Θ F KT = O d d/2 d! = O(2 O(d log d) ). ( 12
)
Note that, in both cases, Θ depends on the ambient dimension but not on δ. Note also that, while similar, the two bounds on the fatness differ by a factor 2 d / √ d. We thus expect that using the Coxeter triangulation as the ambient triangulation will give a smaller output S than the one we obtain using the Freudenthal-Kuhn triangulation. This is confirmed in practice as shown in Section 4 (see Figure 6).

We combine Propositions 22 and 23 in the following theorem. inside σ. To implement the oracle, we first evaluate f at the vertices of σ. We then compute the barycentric coordinates of the (generically unique) point of intersection of the affine hull of σ with H. Lastly, we check whether the barycentric coordinates are all non-negative (to ensure that the intersection point lies inside σ). It follows that the cost of the oracle is the cost of evaluating f at the k + 1 vertices of σ plus the cost of solving a linear system of k equations and k unknowns, which can be done in time O(k 2.375 ).

Dimensionality reduction.

As seen from Proposition 23, the size S of the output of the algorithm, considered as a function of the resolution 1/δ of the triangulation, depends exponentially on m (which is to be expected) and not on d (which is fortunate). Nevertheless, the size of the output depends exponentially on d. This, in particular, means that the sample constructed by the algorithm, although δ-dense, is not guaranteed to be µδ separated for some constant µ. In other words, the output sample is not a net of the manifold.

We can improve on the bound on S by using dimensionality reduction techniques and, specifically, a variant of the celebrated Johnson-Lindenstrauss lemma for manifolds. We depart from our previous worst-case analysis by allowing some approximation factor ε and tolerate a guarantee that holds only with high probability.

Theorem 25 (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 = Ω m ε 2 log 1 ε + 1 ε 2 log Γ δ ,
where Γ is a quantity that depends only on intrinsic properties of M. Let Φ be a random affine subspace of dimension d . Then, with probability

> 1 -η, for all x, y ∈ M (1 -ε) d d ≤ Φx -Φy x -y ≤ (1 + ε) d d .
It follows that the image ΦM of M will be a submanifold of dimension m 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 -m)-simplex σ of the CFK-triangulation of R d , decides whether its inverse image Φ -1 σ,

which is a (d -d )-dimensional flat strip in R d , intersects M or not.
Due to the scaling factor d/d , the resolution of the triangulation in the low dimensional plane 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 (reach and volume) 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. Theorem 24 shows that the S o C G 2 0 2 0 00:

Experimental results

The data structure of Section 2 and the algorithm of Section 3 have been implemented in C ++ . The code is robust and fast, and will be released soon in the GUDHI library [START_REF]GUDHI Editorial Board[END_REF]. Full detail on the implementation, including the implementation of the oracle, will be reported in a companion paper together with experimental results [START_REF] Boissonnat | Tracing isomanifolds of arbitrary dimension and codimension[END_REF]. See also [START_REF] Kachanovich | Manifold meshing using Coxeter triangulations[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.

Data structure.

In Figure 3, we present the time of generating all faces (on the left) and all cofaces (on the right) of various dimensions of simplices in a CFK-triangulation of R 30 using algorithms from Section 2.2. The presented execution time is averaged over 500 tests. Note that both for face and coface generation algorithms, the execution time is proportional to the number of computed elements. On average, these algorithms take time 0.001-0.002 ms per computed face or coface, regardless of the dimensions of the input simplex and of the computed element.

In Figure 4, we further illustrate the particular case of facet and cofacet computation, which is essential in the manifold tracing algorithm. We show the dependency of the execution time on two parameters: the ambient dimension d and the codimension m of the input simplex, which corresponds to the intrinsic dimension of the input manifold in the manifold tracing algorithm.

Manifold tracing algorithm.

We show the performance of our implementation of the manifold tracing algorithm for various ambient and intrinsic dimensions in Figure 5. In Figure 6, we can see that using Coxeter triangulation is beneficial in practice as it produces a smaller output in less time. In

Figure 1 Lemma 4 (Corollary 5 (Corollary 6 .Lemma 7 .

 14567 Figure 1The 2-permutahedron and the ordered partitions associated to its faces.

117 2 . 2 Definition 8 .

 228 Freudhental-Kuhn triangulation118The Freudhental-Kuhn (FK for short) triangulation is obtained from the d-grid, i.e. the 119 unit cubical tessellation of R d that consists of copies of the unit d-cube along the integer 120 lattice Z d . By triangulating each d-cube in the grid in an appropriate way to be described 121 now, we obtain the FK-triangulation of R d . 122 Let x ∈ R d and write z i = x i -x i . We denote by e 1 , . . . , e d the basis 123 vectors and introduce, for reasons that will be clear later, the extra vector e d+1 = -d i=1 e i .

e 1 + e 3 e 3 e 1 + e 2 e 1 + e 3 e 2 e 1 e 2 + e 3 Figure 2

 3121232 Figure 2The permutahedral representation of the simplices in the stars of vertices y and y .

  Given a point x ∈ R d Lemma 9 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.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.

Lemma 16 .

 16 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. 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) from the counting sort.Coface computation. Computing the faces of a simplex σ consisted 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.

Definition 18 .

 18 Let P = {(x i ) ∈ R d+1 | i x i = 0}and consider the d-simplex with vertices u k in P .

Algorithm 1 foreach cofacet φ of τ do 328 foreach facet ρ of φ do 329 if

 1328329 Manifold tracing algorithm 321 input : Triangulation T of R d , manifold M of dimension m, seed point x 0 ∈ M 322 output: Set S of the simplices in T of dimension k = d -m that intersect M 323 Translate T so that x 0 coincides with the barycentre of a k-dimensional face τ 0 in T 324 Initialize the queue Q and the set S with τ 0 325 while the queue Q is not empty do 326 Pop a k-dimensional simplex τ from Q 327 ρ does not lie in S and intersects M then 330 Insert ρ to the queue Q 331

  powerful intersection oracle that not only detects intersections but can also compute intersection points between the simplices in S and M.

Proposition 22 .

 22 The time complexity of the algorithm is O (k2 m I|S|). where I is the time complexity of one call of the intersection oracle and |S| is the size of the output.Proof. The complexity of the initialization is O(d).The complexity of each iteration of the while loop consists of: computing the cofacets of the popped k-dimensional simplex in the queue, computing facets of these cofacets and applying the intersection oracle on each of these facets. From Lemma 7, the number of cofacets is O(2 m ). Each of these cofacets has k + 2 facets. Therefore, for each iteration of the while loop, the algorithm applies the intersection oracle on O(k2 m ) simplices. By using this observation and the complexities in Lemmas 16 and 17, the total time complexity of each iteration of the while loop follows:

  on 1/δ is exponential in m but not in d. Upper bound on |S|. Now, we express |S| in terms of N , d and m. For this, we count the number INC of incidences of the k-dimensional simplices in S and the d-dimensional simplices in N in two ways: τ ∈S |cof(τ, d)| = INC = σ∈N |fac(σ, k) ∩ S|.

Theorem 24 .

 24 The time complexity of the manifold tracing algorithm is 2 O(d log d) I volm(M) δ m , where I is the time complexity of one call of the intersection oracle. Cost of the oracle. The cost of I depends on how the submanifold is given. As an example, consider the case where M is the PL-approximation of the zero set of a function f : R d → R k and assume that evaluating f at any point x ∈ R d can be done in time polynomial in d (which, in particular, is true if each f i is a polynomial in the coordinates of x). Then I depends polynomially on d too. Indeed, consider a k-simplex σ of the triangulation on which we call the intersection oracle and let H denote the m-flat that linearly interpolates f -1 (0)

Figure 3 Theorem 26 ..

 326 Figure 3 On the left: comparison of the execution time of the face and the coface generation algorithm for simplices of various dimensions in a CFK-triangulation of R 30 . Because the average computation time of a face or coface is constant, the presented time is proportional to the number of faces or cofaces of respective simplices.

Figure 4

 4 Figure 4 Execution time of the facet and cofacet computation depending on the dimension d of the triangulation and the codimension m of the input simplex.

Figure 7 , 7 Figure 5

 775 Figure7, we present a PL approximation of a two-dimensional flat torus without boundary embedded in R 10 built by the manifold tracing algorithm. The algorithm can be easily adapted to handle submanifolds with boundary. In Figure8, we present the mesh obtained by our algorithm on a portion of a flat torus embedded in R 4 and cut by a hypersphere.Both surfaces in Figure7and 8 are rotated and translated in their respective ambient spaces for visualization purposes. Note that there is no C 2 embedding of the flat torus in R 3 .

Figure 6 + x 2 2 + x 2 3 - 1 (

 6231 Figure 6 Comparison of the size of the output of the manifold tracing algorithm using two types of the ambient triangulation: 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 1 (x3 + 1) 2 -2x 2 2 = 0 and xi = 0 for i > 3.

S o C G 2 0 2 0 00:4 Sampling and Meshing Submanifolds

  , . . . , ω l ), such that ω i ⊆ T and {ω 1 , . . . , ω l } is a partition of T . The ω i are in a 1 parts, ( a1+1 , . . . , a1+a2 ) is an ordered partition of ω 2 in a 2 parts, . . . , ( a1+...+a k-1 +1 , . . . , a1+...+a k ) is an ordered partition of ω k in a k parts. Note that in Definition 1.3.5, we have c k = |T |. Therefore, for any i 2 {1, . . . , |T |}, there exists an index j 2 {1, . . . , k} such that i 2 {c j + 1, . . . , c j+1 }. Definition 1.3.7 (Refinement). Let ⇡ and ⇢ be two ordered partitions of {1, . . . , d + 1} in k parts and l parts respectively, with k, l 2 {1, . . . , d + 1} such that k  l. We say that ⇢ is a refinement of ⇡ (in l parts), if there exist positive integers a 1 , . . . , a k and ordered partitions ⇢ 1 , . . . , ⇢ k such that ⇢ is a concatenation of ⇢ 1 , . . . , ⇢ k and: In terms of Definition 1.3.7, from ⇢ being the concatenation of ⇢ 1 , . . . , ⇢ k , the number of parts in ⇢ is the sum of the number of parts in ⇢ 1 , . . . , ⇢ k . It follows that the sum of a i for all i 2 {1, . . . , k} is l.

	called the parts. 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]. Definition 3 (Refinement). Let ω and be two ordered partitions of [d + 1] in k parts and l parts respectively, with 1 ≤ k ≤ l ≤ d + 1. We say that is a refinement of ω if there exist positive integers a 1 , . . . , a k such that: ( 1 , . . . , a1 ) is an ordered parti-{1}, {2}, {3} {2}, {1}, {3} {2}, {3}, {1} {3}, {2}, {1} {3}, {1}, {2} {1}, {3}, {2} {1, 2}, {3} {2}, {1, 3} {2, 3}, {1} {3}, {1, 2} {1, 3}, {2} {1}, {2, 3} {1, 2, 3} Figure 1.14: The ordered partitions of {1, 2, 3} that correspond to the faces of a two-dimensional permutahedron. Note that each vertex has three parts, each edge has two parts and the hexagon itself has one part in their corresponding partitions. Remark 1.3.6. ⇢ 1 = ⇢(1), . . . , ⇢(a 1 ) is an ordered partition of ⇡(1) in a 1 parts, ⇢ 2 = ⇢(a 1 + 1), . . . , ⇢(a 1 + a 2 ) is an ordered partition of ⇡(2) in a 2 parts, . . . ⇢ k = ⇢(a 1 + . . . + a k 1 + 1), . . . , ⇢(a 1 + . . . + a k ) is an ordered partition of ⇡(k) in a k parts. We denote the refinement relation by ⇢ v ⇡. Definition 1.3.8 (Set of refinements). Let ⇡ be an ordered partition of {1, . . . , d + 1} in k parts, for some k 2 {1, . . . , d + 1}. Let l 2 {k, . . . , d + 1} be such that k  l. We denote by R(⇡, l) the set of refinements of ⇡ in l parts. tion of ω 1 1.3. The Voronoi diagram of a Coxeter triangulation of type Ãd Remark 1.3.9.

Definition 2 (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 an ordered collection of l subsets ω = (ω 1

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) and the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 754411.

Acknowledgements

We thank Dominique Attali (CNRS), Arijit Ghosh (ISI), Vincent Pilaud (École polytechnique), and Aurélien Alvarez

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.

00:17

A

Proofs for Section 2.4

In this appendix we will prove that the Freudenthal-Kuhn triangulation is the Coxeter triangulations of type Ãd up to a linear transformation. We also prove that the Voronoi cell of a vertex in a Coxeter triangulation of type Ãd is a permutahedron.

For this we need to first recall an equivalent definition of the Coxeter triangulations of type Ãd : Any Coxeter triangulation can be defined as an hyperplane arrangement

where

and R + denotes the set of positive roots of the Coxeter group. We will not recall the general definition of positive roots, which can be found in 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], but use that for Ãd , according to [11, Planche II],

with {s l } the simple roots of the Weyl group A d associated to the triangulation of type Ãd .

For a discussion of the Weyl group we again refer to 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 recall the simple roots of A d , which will be needed in the second proof. The simple roots of A d (in the hyperplane

) are in turn:

see [START_REF] Bourbaki | Lie groups and Lie algebras[END_REF]. We stress that these simple roots can be rescaled, permuted and rotated in the hyperplane P . We note that one can easily rotate R d ⊂ R d+1 given by the first d basis vectors into P . The matrix in SO(d + 1) of this transformation has as first d rows

) and the final row ((1/ √ d + 1) {d+1} ), where c {k} denotes k consecutive coordinates equal to c. We will not use this transformation because it complicates the expressions prohibitively.

In Lemma 9, we have seen that x ∈ R d lies on the face of some simplex with canonical permutahedral representation (ṽ 0 , ω) in the FK-triangulation if and only if either x iṽi 0 =

x jṽj 0 or x iṽi 0 = 0 for some i, j. Note that ṽi 0 , ṽj 0 ∈ Z. Hence we see that

Lemma 27. The Freudenthal-Kuhn triangulation is a hyperplane arrangement

S o C G 2 0 2 0 00:22

Sampling and Meshing Submanifolds

We now define a linear map µ from R d to P by showing how it acts on the basis:

We claim that µ maps E bijectively onto R + . The vector µ(e i ) = r 1,i lies in R + , by construction. For u i,j ∈ E, with i < j, we see that

Hence µ(u i,j ) lies in R + . By reading the previous calculation backwards we see that µ -1 maps each r ∈ R + to a vector in E.

We conclude that µ (bijectively) maps H to H, which completes the proof of Lemma 19.

We now prove the following:

Proposition 28. The Voronoi cell of a Coxeter triangulation of type Ãd is a permutahedron.

Proof. We start by recalling a number of results. In [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] we have seen that the circumcentre of the simplex given in Definition 18 is

) , 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 the simplex of Definition 18 in the hyperplanes of H 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.

The reflection R j,k in a plane that goes through the origin with normal r j,k is given by

We find that

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 [4, 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. Let now 1 ≤ k < j ≤ d be integers such that k + j = d + 1. Then 2 k + 2 j < 2 k-1 + 2 j+1 .

Generalizing this, we see that the sum of the 2 |ωi| -2 is maximal when all |ω i | = 1 except the biggest part which has d -l + 1 elements.