Keywords: Transversalité, Singularités Génériques, Algorithmes Numériques Certifiés, Arithmétique d'Intervalles, Topologie de Courbe Singulière Transversality, Generic Singularities, Certified Numerical Algorithms, Interval Arithmetic, Singular Curve Topology

Entreprendre cette thèse a été une riche expérience, qui m'a révélé non seulement les aspects académiques de ma propre personnalité, mais aussi des caractéristiques intimes, des faiblesses et des forces dont je n'avais jamais pris conscience auparavant. Une thèse de doctorat est souvent décrite comme un travail solitaire; cependant, les noms qui suivent démontrent absolument le contraire. De fait, le chemin vers mon but a été éclairé par l'appui de personnes qui n'ont rien épargné pour m'aider, me motiver et partager des connaissances.

Ce travail a été réalisé avec l'orientation de mes superviseurs: Sylvain Lazard, Guillaume Moroz et Marc Pouget. Avec votre inestimable expertise, la recherche et la méthode allèrent de l'avant. Votre apport perspicace m'incita au raffinement de ma pensée et orienta mon travail vers un niveau supérieur. Merci pour votre amitié, votre patience et votre rigueur académique.

Ma gratitude s'adresse aussi aux membres du jury, qui ont pris le temps de lire cette thèse et ont formulé des commentaires importants et des suggestions, qui m'aident pour l'amélioration de cette contribution et la poussent vers l'avant. Outre l'appui académique, j'ai aussi bénéficié de l'assistance amicale de Joel Veldkamp et de Juan Pedro de Gandt sur le plan linguistique. Vos efforts, bien appréciés, ont enrichi et consolidé la formulation de mon texte, soulignant et renforçant la rigueur et la précision du langage utilisé.

Une profonde gratitude va aussi vers ma famille, pour qui la distance n'a pas été un obstacle aux encouragements et aux preuves d'amour. J'ai aussi bénéficié d'un appui familial de la part de Pierre Kreit et de sa famille, auprès de qui j'ai trouvé de sages conseils et une oreille sympathisante, étant pour moi, en France, comme un second père, en m'offrant, avec Madame Christine Kreit, ce qui me fit sentir chez moi. Le dernier remerciement s'adresse au lecteur de ces pages, quel qu'il soit. Tout commentaire, toute suggestion ou critique seront aussi des contributions pour ce travail. i Résumé L'isolation des points singuliers d'une courbe plane est la première étape vers le calcul de sa topologie. Pour cela, les méthodes numériques sont efficaces mais non certifiées en général. Nous sommes intéressés par le développement d'algorithmes numériques certifiés pour isoler les singularités. Pour ce faire, nous limitons notre attention au cas particulier des courbes planes qui sont des projections de courbes lisses en dimensions superieures. Ce type de courbes apparaît naturellement dans les applications robotiques et la visualisation scientifique. Dans ce cadre, nous montrons que les singularités peuvent être encodées par un système carré et régulier dont les solutions peuvent être isolées avec des méthodes numériques certifiées. Notre analyse est conditionnée par des hypothèses que nous démontrons comme étant génériques en utilisant la théorie de la transversalité; nous fournissons également un semi-algorithme pour vérifier leur validité. Enfin, nous présentons des expériences de visualisation et de robotique, dont certaines ne sont pas accessibles par d'autres méthodes, et discutons de l'efficacité de notre méthode.

This thesis is dedicated to the big family I miss in Syria and to the small one I am making in France.

iii

List of Figures

1.1 Left: Illustration of a configuration space of a 2-degrees-of-freedom robot. In this example, it is a surface embedded in R 3 . Generically, the silhouette of this surface is a (smooth) curve. Right: Projection of the silhouette on the motor space. . . . . . . . . . . . . . . . . 1.2 An approximation of the curve in Figure 1.1 (right) that preserves the topology and the locations of the critical points. The approach used for this goal is to isolate the critical points, and then, to connect them properly using the state-of-the-art methods since the remaining part is smooth. Additional black vertices are added to improve the visualization. The number on each face represents the number of solutions of the direct kinematic problem. . . . . . . . 3.4 Proof of Lemma 4.3.3: The second case where v and v are co-linear.

3.5 Example: The silhouette of the torus in R 3 with the coordinates (x, y, z) (source: [IMP2016a]) . . . . . . . . . . . . . . . . . . . .

3.6

Illustration of Lemma 2.2.17 which is used in the proof of Proposition 3.2.12. Since ϕ is a submersion over Reg(W ) and since j 1 P is, generically, transversal to Reg(W ), the lemma implies that the composition P = ϕ • j 1 P is generically a submersion. . . . . . .

4.6

Comparing the curve C defined in Remark 4.4.5 with C that is defined by x 1 -(x 3 -1) 2 = x 2 -(x 3 -1) 3 = 0 where both of them have the same plane projection which contains an ordinary cusp; however, the multiplicities of S (defined in the same remark) at q differ. In this example, the pair (q, q) satisfies Assumption A - 5 for C , however, it does not for C. . . . . . . . . . . . . . . . . . . . 5.1 An example of Case (ii) for a curve C in R n . The box B contains a smooth part of C that projects to a loop with a node. This implies that two points q 1 and q 2 are in B such that det(M 1 (q 1 )) = 0 and det(M 2 (q 2 )) = 0. Hence, the tangent space plane projection of q 1 (resp. q 2 ) is parallel to x 2 -axis (resp. x 1 -axis). . . . . . . . . . . xiv 5.2 An example of B and B satisfying (2) in Lemma 5.3.1. Their projections have a node in the boundary. The curve C is tangent to the boundaries of both B and B at q and q respectively. The curve and the boxes intersect at q such that q does not have the same plane projection as q and q . . . . . . . . . . . . . . . . . . 8.9 The plane plotting of Isolating_singularities output for S parallel with respect to q 1 and q 2 . . . . . . . . . . . . . . . . . . . . . . . xv 8.10 Plotting the result of applying Isolating_singularities to S serial with respect to x 1 , x 2 , where it plots only red boxes since neither Assumption A 2 nor Assumption A 1 is satisfied. . . . . . . . . . . . 8.11 An example of two different (serial) singular configurations of RRRRR that define the same pose (x 1 , x 2 ). . . . . . . . . . . . . . . . . . 8.12 An illustration of enclosing_singularities output for the serial singularity S serial of the robot RRRRR. The projection is with respect to the (q 1 , q 1 )-plane (the command variables). The certification of the geometric singularities is done only where our strong assumptions are satisfied. The red boxes contain the projections of the (geometric) singular points of S serial . . . . . . . . . . . . . . . . . 9.1 An example where our visualization tool does not provide a correct topology for some precision. A: Our software returns a set of boxes that encloses the projected curve without reporting any node (or unknown-status boxes). B & C: Two different scenarios of the plot in A after enough refinements. . . . . . . . . . . . . . . . . . . .

Introduction

Computational geometry is the branch of computer science that studies theoretical and applied aspects of geometric algorithms. A classical problem in computational geometry is to visualize curves, in particular, plane curves. The importance of curve visualization comes from the fact that curves are very common in engineering, physics and mathematical problems. In general, visualizing a geometric object consists in representing it by an implementable structure which can be handled by computers. In particular, (plane) curve visualization is usually based on approximating a given curve into a piecewise-linear graph that preserves some properties of the curve. A question that might come to mind is: What are the properties that an approximation is supposed to preserve?. The classical properties that are considered are that the approximation should be "close enough" to the original curve and that the original curve and its approximation should have the same topologies, in particular, connected components and self-intersections should not be missed.

Depending on the applications, preserving the topology and some other geometric features may be critical or of interest while in some other applications, this may not be relevant. In this regard, one can find in the literature two main types of approaches for approximating curves:

• Non-certified approaches in the sense that the topology of the approximation does not necessarily coincide with the one of the original curve. In particular, connected components and self-intersection points might be missed or spuriously added. Many numerical non-certified approaches exist in the literature such as marching cubes and continuation methods; see Section 1.3.

• Certified approaches in the sense that they ensure that the topologies of the curve and its approximation are the same. Few certified approaches exist in the literature. Most of these approaches are algebraic, although some numer-ical approaches exist if the input curves satisfy some assumptions. In Section 1.3, we provide more details about these approaches. In particular, we focus on a critical sub-problem that appears in all known approaches of this type, namely, isolating the singularities of plane curves.

In this thesis, we address the problem of certified plane curve approximation. Numerical approaches are usually efficient for approximating curves. However, when the considered curves contain singularities, numerical methods do not, in general, provide certified approximations. Our goal is to design a fast numerical method that is certified for approximating and visualizing plane curves. With this goal in mind, we restrict ourselves to a special case of 2D curves, namely, those that are the projections of smooth curves living in R n , with n > 2. This problem of approximating and visualizing 2D projections of smooth curves in higher dimensions is rather specific and we start by presenting the application in robotics that motivates our work.

Motivation

Before introducing the problem in details, we describe the context that motivates our study. The motivation we present revolves around robotics design and analysis. For this goal, we recall some preliminaries [LP2017]:

An end-effector is the end robotic part that attaches to robotic arms or appendages. The configuration space is the space of possible positions an end-effector may attain. The degrees of freedom of a robot refer to the number of movable joints or links of this robot. Equivalently, the degrees of freedom refer to the dimension of the configuration space. The motor (or command) space is the set of positions that the engines of a robot may attain. The direct kinematic problem refers to determining the position of the end-effector, given the values for the motor variables of a robot. In particular, for a robot with 2 degrees of freedom, the configuration space can be seen as a surface embedded in an n-dimensional space, where n ∈ N represents the number of variables that define and control the robot. Moreover, the engines' possible positions can be represented in R 2 (as a sub-space of R n ).

Figure 1.1 (left) shows how a 2-degrees-of-freedom robot configuration space looks like when n = 3. The first two coordinates x and y represent the engine variables. The third one z represents, for example, the position of an end-effector. This figure is for illustrative purposes and is not meant to refer to any robot. Every point of this surface represents a solution of the direct kinematic problem, given its projection to the (x, y)-plane. The silhouette is the set of points in the surface at which In this example, it is a surface embedded in R 3 . Generically, the silhouette of this surface is a (smooth) curve. Right: Projection of the silhouette on the motor space.

the tangent plane is vertical. Its projection to (x, y)-plane contains the boundary of the engine possible motions. In the neighbourhood of a point on the silhouette, the number of solutions of the direct kinematic problem is likely to change. Figure 1.1 (right) shows the projection of the configuration space's silhouette on the motor space.

Moving the robot from a position to another amounts to convey the engines accordingly. Thus, a path in the motor space is drawn for this motion. If such a path traverses the projection of the silhouette, bad scenarios might occur. For example, suppose that the approximation misses topological information about the silhouette projection such as the part numbered with zero in Figure 1.2 (right). This part surrounds points where the direct kinematic problem has no solution. In this case, paths traversing this part might cause robot damage. One of the approaches used for this problem is based on approximating the silhouette projection into a piecewiselinear graph. This approximation is supposed to preserve several properties, more precisely, the topology and geometric features such as the location of critical points and their types. By critical points we mean those where the tangent line is not welldefined (in this case, the points are called singular) or the tangent line is parallel to one of the axes (smooth critical points). Preserving these properties allows to understand how the silhouette plane projection divides the motor space topologically. This allows to classify the motor space by the number of solutions of the direct kinematic problem. With the help of this information, bad scenarios as the one 1.1 (right) that preserves the topology and the locations of the critical points. The approach used for this goal is to isolate the critical points, and then, to connect them properly using the stateof-the-art methods since the remaining part is smooth. Additional black vertices are added to improve the visualization. The number on each face represents the number of solutions of the direct kinematic problem. mentioned above can be avoided. This type of approximation is also useful when computing paths in the command space for preserving a safe distance between the paths and the boundary of the silhouette. In other words, a topologically-correct approximation of the silhouette projection leads to approximately computing simply connected subsets of the motor space. In every component, the robot is safely movable from and to every two points in that component. More information about this approach is described in Chapter 8.

Recall that, as mentioned in the introduction, if the projections of the silhouette curves contain singularities, numerical methods such as [Sny1992, LOD2001, PV2004, MGGJ2013] fail to ensure a correct topology in general. On the other hand, in the case of algebraic curves, symbolic methods can be used for this goal (see Section 1.3 for more details) but, in general, they are slower compared with numerical methods. However, a specificity of the problem described above is that, although the sil-houette projections generically have singularities, the silhouette curves before the projection are in general smooth. This is why we aim at a certified numerical method for approximating the plane projections of smooth curves in R n .

Although the above example with a configuration space in R 3 is a toy example, other 2-degrees-of freedom robots appear in real-life for different purposes such as plotting, pick-and-place and other industrial tasks. Such robots usually have two motor variables but they tend to have more than one variable for representing the position of an end-effector. Hence, their configuration spaces are surfaces embedded in R n , with n is larger than 3. This explains why we are interested in approximating the projections of silhouette curves in R n and not only in R 3 .

Even with robots of k degrees of freedom, with k > 2, our approach might be helpful by fixing k -2 motors. Then, the motor space is embedded in R 2 and the configuration space is restricted to a surface. We can thus use our approach for approximating the projection of its silhouette. For path planning, this could be used to design safe paths with sequences of motions that only move 2 motors simultaneously.

Problem statement

This section is dedicated to state the problem we consider more formally. We start with explaining the properties that we want to preserve when approximating a possibly-singular plane curve. Then, we describe a natural approach to approximate the targeted type of (plane) curves. Then, we focus on the main sub-problem of this approach, namely, to isolate the singular points of that curve.

Preserved properties when approximating plane curves

Recall that the set of critical points consists of singularities and smooth (nonsingular) points with tangent lines parallel to one of the axis. The curve approximation that we wish is a plane graph where every critical point of the curve is identified by a vertex in that graph. In addition, we previously discussed the properties that we want to certify in order to tackle the robotics problem in Section 1.1.

In the following, we state these properties more formally and justify their significance. To avoid technicality, we describe our purposes by examples. is illustrated in the right-hand side. However, it is clear that this approximation does not look similar to the curve since the relative positions of the connected components are not preserved. This leads us to another interesting property which is beyond topology.

Figure 1.3: Approximating the given curve preserving its topology, but loosing relative positions of the connected components.

To improve the approximation, we require an additional property. This property is about the relative position of the connected components. More formally, we want that the curve and the graph are ambient-isotopic, that is, the curve can be continuously deformed towards the graph such that stopping the deformation at any point gives a curve that has the same topology as the graph.

Geometric properties. By geometric properties we mean the location of critical points. More precisely, for every critical point in the curve and any positive real , we want to find a box in R 2 , with a width at most that contains this point and does not contain any other critical point. This box is called an isolating box. In addition, we also want to determine the type of this point, e.g., smooth critical point, node, cusp, etc. Figure 1.5 shows an approximation of the plane curve mentioned above preserving its geometric properties. Having the critical points isolated, the approximation can be easily improved by adding vertices on critical-point-free parts of the curve. Consequentially, we can control the Hausdorff distance between the curve and its approximation and make it as small as we wish.

From the above discussion, the approach we consider for approximating a plane curve consists in computing an isolating box for every critical point, as well as additional boxes added in the case where there are different edges that have the same endpoints. Then, we identify these boxes as vertices in the approximating Figure 1.4: An approximation of the plane curve that is ambient-isotopic to the given curve. plane graph. This approach allows us to preserve the topological and geometric properties mentioned in the previous paragraph.

Thus, the first step (and the main one) in this approach is to compute isolating boxes for critical points. In the next section we discuss the difficulties of this step.

Isolating singularities

Recall that for implicit plane curves, the critical points are described by systems of equations. In order to detect them efficiently, one may wish to use certified numerical methods such as Newton interval methods [Neu1991]. Recall that these methods require that the system is (a) square, that is it has as many variables as equations, and (b) regular, that is, the system has no multiple solutions; equivalently, the Jacobian matrix at every solution is regular. Generically, critical points can be described by square regular systems in smooth curves. On the other hand, the usual system that encodes singular points (where tangent lines are not welldefined), which is obtained by the Jacobian criteria, does not satisfy theses conditions in general.

For example, take the curve defined by f (x, y) = y 2 -(x 3 + x 2 ) = 0 (Figure 1 that is y 2 -(x 3 + x 2 ) = -(3x 2 + 2x) = 2y = 0 which is not square. Hence, we cannot apply Newton interval methods to certify the correctness of the solutions.

For a system of polynomial equations, symbolic methods are available to certify all solutions of any zero-dimensional system (in particular, to certify curve singularities). However, they show, in general, a slower performance compared with numerical methods; see Section 1.3.2.

From the above discussion, we deduce that the problem of isolating the singularities of a plane curve is itself a challenge towards approximating it. In fact, the main focus of this thesis is to numerically isolate singular points of a special type of plane curves using certified efficient methods.

Formal statement of the problem

In this section, we state the problem more formally and precisely. Before that, we want to determine what we mean by concepts such as "input", "curve" and "solution". In other words, we first mention several related considerations that The first factor to consider is the input of the problem. More precisely, the way in which the curve is given. For example, the curve can be parametrized or, more generally, given by an implicit form, i.e., as the zero locus of a function f (x, y).

In some other problems (as the one mentioned in Section 1.1), the plane curve is given as the projection of a curve in a higher-dimensional space. The plane curves in the previous examples are described by functions. This leads us to another factor, the type of functions used to define that curve, e.g., polynomials, analytic or more general differentiable functions. Moreover, one also considers the domains in which these functions are defined, e.g., real numbers, complex numbers, finite fields, etc. Finally, an important factor is how we want to output singularities (e.g., parametrized, contained in isolating boxes, etc.) and how reliable the output is.

These considerations are parameters of a family of problems, each of which having its own state of the art. In Section 1.3, we discuss the state of the art of the problems that are related to our goal.

Genericity. The formal description of the problem concerning this thesis requires us to mention the concept of genericity. In vague terms, a property is generic if it is satisfied for almost all instances. More formally, the set of C ∞ functions from R n to R n-1 is equipped with the weak topology [Dem2000, §3.9.2] (see Section 2.2.4). The convergence in this topology is defined as follows: For a sequence of C ∞ -functions f i → f , every partial derivative of f i uniformly converges to the corresponding partial derivative of f over every compact subset of R n . A property over C ∞ -functions is generic if it is satisfied by a countable intersection of open dense subsets of functions with respect to this topology.

Formal statement of the problem. Let C be a generic (in particular, smooth; see Lemma 3.2.1) curve defined by a C ∞ -function P from R n to R n-1 , that is, if P maps x = (x 1 , . . . , x n ) → (P 1 (x), . . . , P n-1 (x)), the curve C is the intersection of the n -1 hypersurfaces of the implicit equations P i (x) = 0 for i = 1, . . . , n -1. In particular, C is not necessarily algebraic. Let π C : C → R 2 be the projection sending every point to its first two coordinates. We want to find a set of pairwise disjoint boxes in R 2 (a product of two closed intervals) such that every singular point of π C (C) is contained in a unique box of this set and every box of this set contains exactly one singular point of π C (C).

Plane curves π C (C) of the type described in the latter paragraph appear naturally in visualizing the apparent contour of surfaces which is also related to mechanical design. As mentioned in the motivation (Section 1.1), obtaining information about the silhouette projection to the command space eases the analysis and design of safe motion.

Related work

This Section is dedicated to the state of the art of the methods for approximating plane curves and in particular the plane projections of higher-dimensional curves. We also focus on the sub-problem of detecting their singularities. We are in particular interested in methods that ensure the exactness of the topology and that correctly detect all the singularities.

There are two main approaches for approximating the plane projections of curves in R n . Either the curve is approximated in R n and then projected or it is projected and then approximated.

Most classical numerical methods for approximating curves in R n are based on marching cube (see [GVJ + 2009, Chapter 7]) and continuation methods (see [GVJ + 2009, Chapter 6]). However, even when these approaches may ensure to preserve the topology of the curve in R n [KX1994, FM2007, BL2013, MGGJ2013, TBV2019], they do not ensure to preserve the topology of their projections and we do not detail these methods in this state of the art.

On the other hand, algebraic methods for approximating curves in R n are to the best of our knowledge based on first projecting these curves in the plane (see e.g. [MPS + 2006, Section 3.7]) and they are thus of no interest in the context of our problem.

In summary, there are very few approaches that approximate with a certified topology the plane projections of higher-dimensional curves without first projecting them in the plane. We only know three of such contributions, which we detail in Section 1.3.4: Plantinga and Vegter [PV2006] consider the contour generator projections of implicit surfaces. Delanoue and Lagrange [DL2014] consider the apparent contours of smooth surfaces in R 4 and Imbach et al. [IMP2016b] handle the plane projections of smooth curves in R 3 using a subdivision scheme locally in four dimensions; the latter one is of particular interest for us because our work is a generalization of theirs.

Consider now the approaches that first project the curves in R 2 . When the curves are algebraic, projecting them can be done using elimination theory [CLO2007, Chapter 3]. Then the main problem for approximating the plane curve is to isolate the singularities of the projected curve, which amounts to solving an algebraic zerodimensional system of equations, for which there is a large state of the art (see Section 1.3.2). However, computing the projection (using elimination theory) has a rather high complexity and it is quite time consuming in practice. Furthermore, when the input curves are not algebraic, existing approaches of elimination theory cannot be used.

Our approach is based on the idea of encoding the singularities of the projected curve by a square and regular system that can be solved using certified numerical methods. In fact, the same idea has been already used by Imbach et al. [IMP2016b] for the plane projections of generic curves in R 3 . We provide more details about this work in Section 1.3.4. The same idea was also considered to deal with singular surfaces. More precisely, Diatta et al. [DMP2019] used this idea to encode the curves of singularities that appear on surfaces obtained by projecting generic surfaces in R 4 to R 3 . The rest of this section is organized as follows. In Section 1.3.1 we discuss how to characterize the singularities, given different types of plane curves. Based on that, in Section 1.3.2, we consider the problem of zero-dimensional systems of equations of which singularities detection can be seen as a special case. Then, in Section 1.3.3, we provide a brief state of the art on the isolation of singularities in implicit plane curves. Finally, Section 1.3.4 is dedicated to explain the approaches that are more related to our goal.

Curves and singularities

As mentioned in Section 1.2.3, there are several parameters that must be taken into account when treating the problem of isolating the singular points. When the curve is real or complex and described in the "language of functions" (e.g., as solution set of implicit equations, as a locally or globally parametrized curve, a projection of a curve embedded in a higher dimension, etc.), the Jacobian criteria can be used to encode singularities as the solution set of a zero-dimensional system of equations. If the curve is parametrized by (f 1 (t), f 2 (t)), then every singular point corresponds to a t-value with (f 1 (t), f 2 (t)) = (0, 0). If the curve is given by f (x, y) = 0, with f is differentiable, then the singularities are the solutions set of the system f (x, y) = ∂f ∂x (x, y) = ∂f ∂y (x, y) = 0. For the case where the curve is the projection of an implicit curve in a higher dimension, it is more complicated. As mentioned in the preamble of this section, one may project the curve to R 2 to compute an implicit equation that defines the projection by which the singularities can be characterized. Another way to characterize the singularities can be done in terms of points in the upper curve as in [IMP2016b, DL2014] which involves zero-dimensional systems of equations.

From the foregoing, our considered problem can be seen as a special case of solving a zero-dimensional system of equations. This general problem is fundamental in algebra and has itself parameters that change many aspects. So, in addition to the state of the art of singularities detection, it is useful to go through the history of solving a zero-dimensional system of equations. This allows us to better understand the different available approaches of singularities detection. In general, two main types of approaches are used to deal with both problems. In the following section, we clarify these types. We discuss their strategies, pros and cons that interest us.

Zero-dimensional systems

By a zero-dimensional system we mean a system of which the solutions set forms a zero-dimensional manifold. This is a very active area of research involving, in general, two main types of approaches. These approaches are numerical and symbolic. In the following, we explain these types in details.

Numerical methods

Numerical methods are tools to solve mathematical problems where the output is an approximated form of the targeted object. In the context of solving system of equations (in particular, isolating singularities of plane curves), the basic strategy in numerical methods is to construct sequences of objects that converge towards the solutions. This gives rise to an important question, that is whether the convergence is always guaranteed. Unfortunately, the answer is no in general. Moreover, the solutions are usually sought in a given area of interest. This makes searching for roots in precise regions efficient. However, if the goal is finding all solutions, numerical methods do not give a complete answer, unless bounds on the solutions are given. On the other hand, numerical methods provide fairly fast approaches compared with the other category (which we introduce next). Moreover, in terms of the type of equations, numerical methods can solve a wide range of systems such as analytic and differentiable.

One can find in the literature classical approaches based on Newton (-Raphson) or Fixed-Point methods; see for example [Epp2013, ( §) 7.8] for details. These methods are not certified in general. However, under some conditions, the convergence can be guaranteed. Concerning our goals in this thesis, we focus on certified methods. An important type of numerical methods is subdivision-based methods; see for example [EH2004, SP1993, GS2001b] and the references therein. Subdivision can be seen as a generalization of binary search. Starting from a given domain, the subdivision consists in recursively bisecting the area until existence and uniqueness criteria are satisfied or the non-existence of solutions is guaranteed. One can distinguish two main types of subdivision: algebraic and box-function-based. Algebraic subdivision handles systems of polynomials. One exploits representations of polynomials such as Bernstein form or B-splines [ZSP1993, GS2001b, GS2001a, MP2009]. On the other hand, boxfunction-based subdivision [Neu1991, PVHK1997, Kea2013] is based on interval evaluations of the given functions (and possibly their derivatives) which supports a broader class of functions.

In these approaches, different tests are used. In the case of differentiable functions of which box functions are given (i.e., interval enclosures of the function evaluations; see Definition 5.1.1), some approaches are based on interval computation techniques such as interval evaluation, Newton interval method, Krawzcyk operator [Neu1991, MKC2009]. These operators offer existence, uniqueness and exclusion tests (see [Neu1991, Theorems 5.1.7 & 5.1.8]). Other approaches use tests such as fixed point theorem of Brouwer, in the form of Miranda's Theorem [Kul1997]. For example, in [XY2019], Miranda and Jacobian test are used to verify the existence and uniqueness respectively. In addition, the exclusion test of the previous paper is based on interval evaluation.

Another type of numerical methods is based on homotopy continuation; see for example [BL2013,VDH2015,TBV2019]. Using a system S with known solutions, homotopy-based methods build a continuous path between the solutions of a given system S and the solutions of S .

Symbolic methods

Roughly speaking, symbolic methods are algorithms with exact input and output based on a formal language and algebraic-combinatorial approaches. In the case of solving a zero-dimensional system of equations (in particular, isolating singularities), these methods are usually based on elimination strategies. Thus, they can be seen as generalizations of the classical Gauss elimination method. The basic idea of this type of methods is, first, to transfer the given problem to a lower dimensional one using algebraic-combinatorial methods, returning a formal representation of the solutions. This step is usually called a symbolic or a projection step. The second step is to transfer this characterization of solutions into an approximated numerical presentation. This step is called numerical or lifting step. In contrast to general numerical methods, symbolic methods are certified and provide a complete answer about the solution set. Moreover, they do not require initial conditions to guarantee the correctness or completeness of the output. On the other hand, in most well-conditioned problems, the techniques used in the projection step are usually slow. compared with the numerical ones. This makes symbolic methods unpractical for systems of high dimensions and degrees. Moreover, most algebraic techniques apply only for polynomials which limits the generality of such approaches.

In general, two main types of symbolic methods are used for zero-dimensional systems: those that are based on resultant techniques or Gröbner basis. Using these approaches, solving a zero-dimensional system amounts to, roughly speaking, computing a formal representation of the solutions such as a univariate representation or a triangular representation. To stay in the context of our thesis, we present in the following a short survey of the contributions considering mainly the bivariate case. For the general zero-dimensional case, we refer the reader to [CLO2007] and [RZ2009].

Resultant techniques. The resultant of two univariate polynomials is a polynomial expression of their coefficients that is zero if and only if the two polynomials have a common root. This idea can be generalized to higher dimensions. In the latter case, the resultant is used for computing the projection of varieties. For solving zero-dimensional systems (in particular isolating singularities), the usual approach using resultant is to compute the projection of solutions along several directions. Then, to obtain candidates which are then filtered until we get only the solutions. In the context of implicit plane curves, methods based on genericity are presented in [GVK1996,SW2005,CGL2009,Vil2018]. Other papers use resultant techniques without generic position assumption such as [DET2009, BES2011].

Gröbner basis. The Gröbner basis technique was introduced by Bruno Buchberger [Buc1965].1 Given a system of polynomial equations, a Gröbner basis of this system is a set of polynomials satisfying the following: A polynomial f vanishes in all solutions of the system iff it is divisible by the Gröbner basis, that is, the reminder of f , over the Gröbner basis polynomials, by the Euclidean division, is zero. In additional to zero-dimensional systems, Gröbner basis techniques are used in many applications in computer algebra; see for example [GM1989, Laz1992, Sva2014] and the references therein.

Based on the previous techniques, the following representations can be used to solve zero-dimensional systems:

Univariate representation. Given a zero-dimensional system of polynomial equations, this representation consists of, first, a univariate polynomial. The roots of this polynomial are in one-to-one correspondence with the given system solutions. Second, a set of univariate rational mappings that parameterize the coordinates of the solutions, that is, the mappings send the roots of the univariate polynomials to the solutions of the system. This type of representation started in end of 19th century with Kronecker's work. This approach usually consists of two steps: The first one is to compute a polynomial (called separating polynomial) that takes pairwise different values when evaluated in distinct solutions. Second, to compute the polynomials defining the univariate representation.

Recently, this type of techniques has been widely studied for both solving zerodimensional polynomial systems and singularity detection of implicit plane curves; see for example [ABRW1996, Rou1998, BSS2001, Bou2014, BLM + 2016] and the references therein.

Triangular representation. The basic idea of this method is to recursively rewrite the system that we want to solve in a triangular way, that is, in an ordered list of equations such that every equation has one variable less than the one before. Hence, in the manipulated system, there is a univariate equation which can be solved numerically. This allows to analogously compute approximations of the other variables. In the bivariate case, this amounts to saying that the manipulated system is of the shape f 1 (x, y) = f 2 (x) = 0. The function f 2 is solved using univariate numerical methods [Epp2013, §3]. Then, f 1 is solved in the same way using the computed approximations of f 2 -solutions; see for example [Laz1992, Kal1993, WCB2001, vdHL2018] and the references therein for more details.

Singularities of implicit plane curves

By an implicit plane curve we mean the zero locus of a bivariate function f (x, y). To approximate such curves, the literature is rich of uncertified numerical methods that compute approximations of given curves without any guarantee in the topology. Since this type of methods is out of our interests for this thesis, we simply refer the reader to [GVJ + 2009] for more details.

For approximating plane curves in a certified way, one can find in the literature many approaches that are based on subdivision such as in [Sny1992] or homotopy such as in [BL2013]. These approaches are certified under the smoothness condition. Some software is also available to certify smooth curves such as Ibexsolve [Nin2015]. This software provides a tool to enclose smooth plane curves, even for some non-algebraic curves, more precisely, for those that are defined by the Minibex language (e.g., sin, cos, exp. . . ).2 For singular curves, ISO-TOP [CLP + 2010] can be used to certify the topology of algebraic curves.

Concerning the problem of isolating singularities, recall that, for plane curves, the singularities are encoded by the system f (x, y) = ∂f ∂x (x, y) = ∂f ∂y (x, y) = 0. This system does not satisfy the conditions required for applying certified numerical methods such as the interval Newton methods.

A numerical approach to solve this issue is to use combinations of these equations to translate this system into a square one.

Nonregular solutions can be handled through deflation systems (see for example [OWM1983, LVZ2006, MM2011, HMS2015]). However, spurious solutions may appear in this process. We are not aware about any approach that finds the singularities of a general implicit plane curve numerically in a certified way. Recall that, in Section 1.3.2, we discussed the symbolic methods for isolating the singularities of plane curves in the context of solving zero-dimensional systems of equations.

Singularities of the plane projections of implicit space curves

As mentioned in Section 1.1, such curves appear naturally in real-life problems as in robotics. In the polynomial case, for curves given as the intersection of two surfaces in R 3 , the plane projection can be defined by the resultant of the polynomials defining the surfaces. Its singular locus can be related to the first sub-resultant. Imbach et al. [IMP2016a] use the latter facts and present a deflation-based approach to generically isolate the singularities of the plane projection (i.e., the zero set of the resultant). This approach involves obtaining a square regular system that encodes the singularities. The idea of this system is based on a relation between the singularities of the plane projection and sub-resultants under some generic assumptions. In other words, under some generic assumptions presented in [IMP2016a, §2], the authors provide a square regular system [IMP2016a, Theorem 2], based on subresultants that can be used to numerically isolate singularities in a certified way.

A different approach is presented in [IMP2016b] for more general curves, namely, those that are defined by analytic functions. Similarly to [IMP2016a], the goal of the latter work is to numerically isolate the singularities of the plane projection of a curve in R 3 . On the other hand, the approach in [IMP2016b] uses information about the projection function Π xy from the space curve onto R 2 . More precisely, under some generic assumptions [IMP2016b, §2.1], the authors show that a node (resp. an ordinary cusp) appears in the plane projection exactly where Π xy is not injective (resp. not immersion; see Definition 2.2.6). Based on these statements, a regular square system that encodes singularities is presented. This system is called the Ball system and is solved using subdivision approaches. Comparing with the system presented on [IMP2016a], the Ball system is of higher dimension which may mean that it is more costly to solve. On the other hand, it is more general since it holds for analytic curves. An algorithmic version of this approach in addition to an implementation (for polynomials with interval coefficients) are presented in [IMP2018]. Even though subdivision approaches may suffer in practice from the curse of dimensionality, Imbach et al. observe experimentally that, for algebraic curves, their approach is more efficient than computing the implicit equation of the projected plane curves and its singularities using symbolic methods. In fact, our work in this thesis can be seen as a generalization of this approach for generic space curves in R n , for n > 2 defined by C ∞ functions.

Plantinga and Vegter [PV2006] provided for the first time explicit conditions on implicit surfaces in R 3 to exhibit generic singularities in the contour generators (the curves on the surfaces which separate front-facing from back-facing regions) and their projections. They presented a certified algorithm to approximate the contour generator and its projection. The algorithm is based on a theoretical analysis which provides local models of the singularities. Delanoue and Lagrange [DL2014], compute the topology of the apparent contour (the set of critical values) of a smooth mapping from R 2 to R 2 . The reader can notice the connection of this work with our considered problem by the following: The graph of a generic function f : (x, y) → (X = f 1 (x, y), Y = f 2 (x, y)) can be seen as a parametric surface in R 4 , defined by {(x, y, X, Y ) ∈ R 4 | X = f 1 (x, y), Y = f 2 (x, y)}. The silhouette is the set of those points in the surface where the projection of the tangent onto R 2 is not surjective. The plane projection of the silhouette is equal to the apparent contour. This approach is based on result proved by Whitney that, generically, critical points of a smooth map are folds or cusps [Whi1955]. Using interval Newton method, Delanoue and Lagrange solve two systems of equations which represent ordinary cusp and nodes respectively. The first one is based on the idea that f restricted to the apparent contour around an ordinary cusp is equal (up to an equivalence) to (x, y) → (x, xy + y 3 ). The second system is based on characterizing distinct points of the domain that have the same image under f .

Contribution

In this section we give a detailed overview of the results we obtained. Recall that we are interested in the problem of isolating the singularities of the plane projection of a generic curve C of R n , with n > 2. Our approach is based on the following: First, we find proper conditions on C that lead to characterizing the singularity in a useful way. These conditions are general enough in the sense that we prove that they are generic. Hence, after finding these conditions, we prove their genericity. Then, under these conditions, we introduce a regular-square system that encodes the singularities of π C (C). Based on this system, we develop a (semi) algorithm that, first, checks whether a given curve C satisfies our assumptions. Second, it returns a set (or a superset) of boxes in R 2 that isolates the singularities of π C (C). In addition, we provide an implementation of our (semi) algorithm and we examine its performance using examples of implicit curves and an application in robotics.

The main novelty of our work is that we consider types of curves for which no certified methods used to exist. Furthermore, our approach allows us to handle classes of curves that were in practice not reachable by other methods in a certified way. Namely, we consider curves defined by C ∞ -functions and not only algebraic or even analytic functions. Moreover, for curves that live in R n for n larger than 3, no certified numerical methods existed even for algebraic curves. For such types of algebraic curves, our preliminary experiments show that our approach allows us to handle, in practice, curves that are unreachable by classical algebraic certified approaches.

We give in Section 1.4.1 a summary of the assumptions we consider. Then, in Section 1.4.2 we explain our idea for encoding the singularities in a square and regular system, which we call the Ball system. Section 1.4.3 is dedicated to briefly describing the specifications of the algorithm we designed for checking the assumptions and isolating the singularities of the projected curve. In addition, we show how we overcame the curse of dimensionality when solving the Ball system. Finally in Section 1.4.4, we outline the properties of our implementation and the experiments we perform to measure its efficiency.

Generic assumptions

Our goal is to design a numerical algorithm for isolating the singularities that appear in the plane projection of a curve C in R n . Numerical algorithms usually cannot handle degenerate cases, that is, singularities in our context. However, under some assumptions on C (which are formally defined in Section 3.1.2), we succeed to handle in a certified way some, although not all, singularities of the projection. Namely, we require that the singularities are "generic", that is, only nodes can appear in the projection. Our other assumptions on C are, roughly speaking, that it is smooth, that its projection only has a discrete set of singularities, and that at most two points of C project on each singularity. We call this set of assumptions strong. We prove in Section 3.2 that these assumptions are generically satisfied.

In addition, we consider a weaker assumption that relaxes the constraint that only nodes can appear. This assumption allows for ordinary cusps in π C (C). Our motivation for also considering this weak assumption is dual. First, our certified algorithm for isolating the singularities of the projection of curves satisfying the strong assumptions also works, to some extent, if only the weak assumptions hold: namely, it outputs a superset of the isolating boxes of the singularities. Second, we conjecture that our weak assumptions are satisfied by silhouette curves of generic surfaces (see Proposition 3.2.17 and Conjecture 3.2.18). A detailed explanation about the assumptions and their genericity are presented in Chapter 3.

Enclosing singularities by the Ball system

We consider here our set of weak assumptions, that is, nodes and cusps are allowed in the projected curve π C (C). Our goal is to encode the singularities of π C (C) by a square and regular (see Definition 4.4.1) system so that it is solvable with certified numerical methods.

Intuitively, a node appears when two points of C project to the same point and a cusp appears when projecting a point with a tangent line orthogonal to the projection plane (see Figure 1.7). The idea to encode the nodes is to design a system whose variables are the coordinates of two different points in R n constrained to be on C and so that they have the same plane projection. To encode a cusp, we design a system whose variables are the coordinates of one point in R n constrained to be on C with a tangent orthogonal to the projection plane. Furthermore, we need systems that are square and regular for using certified numerical solvers. To solve this issue and to gather the two systems into a single one, we first parameterize two different points of C with the same projection by (x 1 , x 2 , y + r √ t) and (x 1 , x 2 , y -r √ t), with x 1 , x 2 , t ∈ R, y, r in R n-2 and r = 1, where r denotes the Euclidean norm of r. Hence, the Ball system is of size 2n -1. Then, given any function f from R n to R so that f = 0 is one of the n -1 hypersurfaces that define C, we introduce (Definition 4.1.1) two smooth functions S • f and D • f . When t > 0, they return, roughly speaking, the arithmetic mean and difference of f at the above two points, hence they both vanish if and only if the two points are on the hypersurface f = 0. When t = 0, the two points coincide and S • f and D • f return, roughly speaking, f evaluated at this point and the gradient of f (at that point) scalar the "vertical" vector (0, 0, r); hence, they both vanish if and only if the point is on the hypersurface f = 0 and its tangent hyperplane is normal to the plane of projection. It follows that given a curve defined by P 1 = • • • = P n-1 = 0, the solutions of the so-called Ball system of all S • P i = D • P i = 0 is the set of points on the curve that project to nodes and cusps (Theorem 4.1.3).

Note that we consider √ t instead of t in the parameterization (x 1 , x 2 , y ± r √ t) for ensuring the regularity of the Ball system when t = 0 (because this ensures that the linear term of the Taylor expansion of D • f , with respect to t, does not vanish). We refer the reader to Chapter 4 for a more formal explanation.

Algorithmic contribution

Based on interval arithmetic [MKC2009], we provide a semi-algorithm to check whether a given implicit curve in R n satisfies the weak (or strong) assumptions. More precisely, the semi-algorithm terminates if and only if the assumptions are satisfied. Moreover, under the strong assumptions the semi-algorithm provides a set of pairwise disjoint boxes in R 2 each of which contains exactly one singular point (in this case, it is a node) of π C (C) such that the boxes union contains all singularities. On the other hand, under the weak assumptions, the semi-algorithm returns a superset of the strong-assumption-case output. A full description of the semi-algorithm we designed is in Chapter 5.

Pruning the domain of the Ball system. As mentioned before, the Ball system is of size 2n -1. This is because the Ball system encodes (a) two distinct points in the case where the singular point is a node (see Figure 1.8, left) or (b) a point and a tangent vector for the case of an ordinary cusp. It is arguable whether solving a system of size 2n -1 for detecting singularities of a 2D curve is a good idea

C is smooth π C (C) Cusps Node (x 1 , x 2 , y) (0, 0, r √ t) y ∈ R n-2 x 1 x 2 (0, 0, -r √ t)
(0, 0, r) from an efficiency point of view. In fact, we succeed to overcome the curse of dimensionality by pruning the domain on which we solve the Ball system.

The idea is that we compute a set of n-dimensional boxes that cover C. Then, every singular point of π C (C) is contained, roughly speaking, in either (a) the plane projection of a box where the curve has (possibly equal) x 1 -critical points and x 2critical points or (b) a non-empty plane projection of two disjoint boxes. Then, we map these boxes to boxes in R 2n-1 ,3 which are the restricted domain on which we solve the Ball system (see Figure 1.8, right). A detailed description of the pruning in addition to other improvements are described in Section 5.3.

Implementation, experiments and applications

We provide a software Isolating_singularities4 written in Python that is an implementation of our semi-algorithm mentioned in Section 1.4.1. This software uses black boxes such as Python-FLINT [Joh2012b] and Ibexsolve [Nin2015] to deal with interval-arithmetic objects and to solve regular systems of non-linear equations. In addition, our software provides a way to visualize the projection of a space curve into R 2 or R 3 . We present our software in Chapter 6 and experiments and discussions on its performance in Chapter 7. We terminate the summary of our contribution with exhibiting results regarding an application in robotics. Within this frame, we study the plane projection of the parallel and serial singularities (not to be confused with the geometric singularity; see Definition 8.1.1 (f)) of a parallel 2-degrees-of-freedom robot. In the considered example, the parallel and serial singularities are curves in R 4 defined by analytic functions. With respect to different plane projections, we compute certified isolating boxes of singular (the geometric sense) points in the plane projection of the parallel and serial singularities (as long as our assumptions are satisfied). In the end, we provide 2D and 3D plots to visualize the results. See Chapter 8 for more details.

Part I

Singularity modelling

This part presents the main content of this thesis. It contains the theoretical aspects of our contribution and the results that prove the correctness of our approach. This part consists of local and global analysis of the plane projection curve in addition to (semi-) algorithms that are based on our theoretical results.

Before starting the analysis, we recall some facts in algebra and geometry in Chapter 2. Then, in Chapter 3, we introduce a set of conditions on curves in R n that guarantee the correctness of our approach. We prove the genericity of these conditions discussing, in the end of the chapter, the same conditions on the silhouette of a generic surface. In Chapter 4, we introduce the Ball system that encodes the singularities the plane projection of a space curve satisfying our conditions. We show, in addition, the relation between the regularity (Definition 4.4.1) of the Ball system and the generic conditions.

The algorithmic part consists of Chapter 5, where we provide the semi-algorithms for isolating the singularities of the plane projection of the space curve. We start the chapter with preliminaries in interval analysis. After presenting the semialgorithms, we prove in Section 5.2 the correctness in addition to determining the necessary and sufficient conditions for termination. In Section 5.3, we explain improvements to the semi-algorithms preparing them to the implementation stage. This chapter is dedicated to recall some concepts and results in algebra and geometry. Most of the material in this chapter is well-known and can be found in any introductory textbook on commutative algebra, algebraic and differential geometry. We suggest to the reader the references [CLO2007, CLO2005, Dem2000]. The rest such as transversality theorems are less known.

We start this chapter with common properties of the ring of real polynomial and determinantal varieties. Then, we recall some results about manifolds and singularity theory. After that, we exhibit the concept of the multiplicity of non-linear system of equations. This concept is a generalization of the intersection multiplicity of the algebraic case. We terminate the chapter with transversality theorems which have a direct connection of our assumptions.

Real algebraic geometry

In this section, we recall concepts in algebraic geometry that we use in different places in this thesis. We assume that the reader is familiar with elementary notions and well-known results of commutative algebra such as polynomial rings, localization, ideals, Gröbner Bases, etc, and algebraic geometry such as affine varieties, Hilbert's Nullstellensatz, etc. We suggest to the reader the reference [CLO2007], in particular, Chapters 1, 2, 4, 5 and 9.

The ring of real polynomials

For a positive integer n and x = (x 1 , . . . , x n ), we consider the polynomial ring R[x] over the field of real numbers R. In the following, we present the definition of the multiplicity of a zero-dimensional ideal over R[x] at a solution. This notion will be used to count the preimages of points in R 2 under π C (see Section 1.4.1). Loosely speaking, the multiplicity m at a solution a in this setup is a generalization of the univariate case where a is a repeated solution m times. For a point q ∈ R n , consider the ring R[x] q that is the localization of R[x] at q. Definition 2.1.1 ([CLO2005, Definition 4.2.1]). For integers m n 1, let G = (g 1 (x), . . . , g m (x)) be a polynomial function from R n to R m and q be a solution of the system {G = 0}. Define the intersection multiplicity of q in the system {G = 0} (or equivalently the multiplicity of the system {G = 0} at q) to be the dimension of the real vector space R[x]q I G , where I G is the ideal generated by the set

{ g 1 1 , . . . , gm 1 } in R[x] q .
The previous definition is classical for the algebraic case. However, we are interested in curves defined as the zero locus of smooth functions. For this goal, we consider a more general definition in 2.2.10 for a system S = {f 1 (x) = • • • = f m (x) = 0} with f i is differentiable infinitely many times which is, unsurprisingly, more complicated and its computability is less clear. Definition 2.1.1 will be thus used when it is equivalent to Definition 2.2.10.

In some places of this thesis where we compute the multiplicity, we reduce the problem of computing the dimension of R[x]q I G to computing the dimension of another quotient that is easer to handle. The following theorem is in this regards: Theorem 2.1.2 ([CLO2005, Theorem 4.4.3, P. 177] ). Consider the real polynomial ring R[x] and its localization R[x] q at some point q ∈ R n . Consider a local monomial ordering > over R[x] q . Let I be an ideal in R[x] q and LT (I) be the ideal generated by the leading monomials of I. Then, the following statements are equivalent:

(a) The dimension of the R-vector R[x]q I space is finite.

(b) The dimension of the R-vector R[x]q LT (I) space is finite.

Furthermore, when any of these conditions is satisfied, we have:

(c) The equality dim( R[x]q I ) = dim( R[x]q LT (I) ) holds.

Determinantal varieties

We study in this section a large and important class of varieties, those whose equations take the form of the minors of a matrix. This concept is useful in Section 3.2.2 to prove the genericity of properties for a silhouette curve, where those properties are based on determinantal relations.

Definition 2.1.3. Let M (m, n) be the vector space of real matrices of size m × n and r be a positive integer such that r < min{n, m}. The determinantal variety, M r (m, n), is the set of matrices in M (m, n) that have rank less than r + 1.

Proposition 2.1.4 ([BV1988, Proposition 1.A.1.1]). Let M r (m, n) be a determi- nantal variety in M (m, n).
The following statements hold:

(a) M r (m, n) is an irreducible variety in M (m, n). (b) M r (m, n) is of dimension r(n + m -r). (c) The singular locus of M r (m, n) is M r-1 (m, n).
Lemma 2.1.5 ([Bôc1964, §XIV.61 Theorem 1]). Let n 2 be an integer, {x ij } 1 j,i n be a set of n 2 variables and C[x ij ] 1 j,i n be the ring of complex polynomials with the variables {x ij }. Then, the determinant of the matrix

(x ij ) 1 i,j n is an irreducible polynomial in C[x ij ] 1 j,i n .

Differential geometry

In this section, we study real C ∞ (or smooth) functions, i.e., those that are differentiable infinity many times. Since the main goal of our work is to deal with the singular points of plane curves given in terms of smooth functions, it makes sense to study the behavior of this type of functions (and their zero locus) locally.

Smooth functions and manifolds

In this section we state the properties we need of the type of functions that interest us, namely, C ∞ functions. More formally, let n and m be positive integers. A function f defined on an open subset of R n to R m is called C ∞ (or smooth) at a if it is differentiable infinitely many times at a, that is, all partial derivatives of f at a of all orders exist. For an open subset U ⊆ R n , the set of smooth functions from U to R m is denoted by C ∞ (U, R m ). In this section, we state useful statements about C ∞ functions and their zero locus.

Definitions 2.2.1 (Regular and singular points [Dem2000, Definition 2.2.2]). Let m 1 be an integer, V be a subset of R m and p ∈ V . We call p a regular (or smooth) point of

V if V is a sub-manifold at p, that is, there exist a neighbourhood W of p in R m , an integer k > 0 and k smooth functions ϕ 1 , . . . , ϕ k defined over W , such that V ∩ W is the set of solutions of {ϕ 1 (x) = • • • = ϕ k (x) = 0} in W
and the rank of the matrix

   ∂ϕ 1 ∂x 1 . . . ∂ϕ 1 ∂xm . . . . . . . . . ∂ϕ k ∂x 1 . . . ∂ϕ k ∂xm  
 , evaluated at q, is k. We call this matrix the Jacobian matrix of the system {ϕ 1

(x) = • • • = ϕ k (x) = 0}
and we denote it by J (ϕ 1 ,...,ϕ k ) . If q is not a regular point of V , we call it a singular point (or singularity). If all points in V are regular, then V is called regular or smooth. Otherwise, V is called singular.

For ϕ = (ϕ 1 , . . . , ϕ k ) ∈ C ∞ (R n , R k ),
we denote by T q ϕ its derivative (also known as the tangent map) at the point q. Note that the Jacobian matrix J ϕ = J (ϕ 1 ,...,ϕ k ) is the expression of the derivative in the canonical bases of R n and R k .

A function f : X → Y between two manifolds is called diffeomorphism if f is a bijection and T q f is full rank at every q in X. Two manifolds are called diffeomorphic if there exist a diffeomorphism from one of them to the other.

The following Theorem shows the compatibility between the dimension as a (smooth) real variety and the dimension as a manifold:

Theorem 2.2.1. [BCR1998, (Proposition 2.8.14)] Let V ⊆ R n be
a manifold that is also an affine variety. Then its dimensions as a manifold and as a variety coincide.

The following theorem by Whitney represents a nice property of even smooth functions. We use this theorem to rewrite functions locally and parametrize them in terms of variables of an even power.

Theorem 2.2.2 ([Whi1943, Theorem 1 & 2]). Let f be an even (resp. odd) smooth function, then there exists a smooth function g such that f

(x) = g(x 2 ) (resp. f (x) = x • g(x 2 )). U W f (x, y) = 0 g(x, y) = 0 q Figure 2.1: An example in R 2
where f and g define the same germ.

Germs of smooth functions

As mentioned before, we are interested in the local behaviour of a smooth function that are defined in R n , that is, we want to study a function that is defined in a neighbourhood of a point q ∈ R n . For this reason, we can assume in this section, without loss of generality, that q = 0 ∈ R n . Definition 2.2.3. [Dem2000, §5.2.1] Let U be an open sub set of R n with q ∈ U and f, g : U → R be two smooth functions. The functions f, g define the same germ at q, if there exists an open subset W contained in U with q ∈ W such that f W = g W , that is, f and g restricted to W coincide. Consequentially, the zero locus of f, g coincides in W (see Figure 2.1) Definition 2.2.4 ([AGZV2012, §17.1]). For i ∈ {1, 2}, let C i be a plane curve defined in a neighbourhood U i ⊂ R 2 of p i by the zero set of a smooth function f i . The pairs (p 1 , C 1 ) and (p 2 , C 2 ) are equivalent, and thus define the same plane curve singularity, if there exists a diffeomorphism ϕ from

U 1 to U 2 such that f 1 = f 2 • ϕ and ϕ(p 1 ) = p 2 .
In particular, a singularity is of type A k if the curve is locally defined at the origin by the zero set of the function x 2 -y k+1 . As important special cases, A 1 is called a node singularity and A 2 is called an ordinary cusp singularity, see Figure 2.2. We are going to assume later that the plane curve that we consider has only nodes or ordinary cusps as singularities. After that, we prove that this assumption is "general enough", more formally, generic; see Chapter 3. Remark 2.2.5. It is worthy to notice that a curve C has an ordinary cusp at a point p if C can be locally parametrized with (z 2 , z 3 ) and p corresponds to the value z = 0. This remark is helpful to characterize ordinary cusps in Chapter 4.

Definition 2.2.6. [GG1973, I.2.3] Let X and Y be differentiable manifolds. Let f : X → Y be a differentiable function. Suppose that T q f is full rank. Then,

• if dim(X) dim(Y ), f is an immersion at q, • if dim(X) dim(Y ) and f (X) is diffeomorphic to X, f is embedding • if dim(X) dim(Y ), f is an submersion at q,
• if for every point q ∈ X, f is immersion (resp. submersion) at q, then f is an immersion (submersion).

Proposition 2.2.7. [Dem2000, Proposition 2.9.6] Preserving the same notation as in Definition 2.2.6, if f is an immersion at q ∈ X, there exists an open subset U in X containing q such that the restriction of f to U is embedding.

The multiplicity of smooth systems of equations

We defined earlier the multiplicity of zero-dimensional systems of polynomials.

Our goal now is to define the same notion for smooth functions. We first recall the definition of multiplicity in the univariate case before generalizing it to higher dimensions.

Definition 2.2.8. Let f be a real smooth function at a ∈ R. The order of f at a is the integer

ord a (f (x)) = min{k ∈ N | ∂ k f ∂x k (a) = 0} if it exists, otherwise ord a (f (x)) = ∞.
For the case a = 0, we write for simplicity ord(f ) = ord a (f ).

The following definition is introduced in [DLZ2011]. We are going to use the same notion of the previous paper:

2.2. Differential geometry Notation 2.2.9. Consider x = (x 1 , . . . , x n ), for k = (k 1 , . . . , k n ) ∈ Z n and q ∈ R n , define the differential operator ∂ k[q] : C ∞ (U, R) → R that sends f to 1 k 1 !...kn! • ∂ k 1 +•••+kn f ∂x k 1 1 ...∂x kn n (q) if k 0 ∈ Z n
and 0 otherwise. From now on, a linear combination of differential operator is called differential functional. The set of all differential functionals forms a vector space denoted by D q over R. Define the linear antiderivative operator φ j that maps

∂ k ∂x k 1 1 ...∂x k j j ...∂x kn n to ∂ k-1 ∂x k 1 1 ...∂x k j -1 j ...∂x hn n if k j > 0 and to ∂ 0 x 0 j otherwise, where k = n i=1 k i . Let S = {f 1 (x) = • • • = f m (x)
= 0} be a system of smooth equations, q be a solution of S and k be a non-negative integer, we define the dual space of rank k, denoted by D k q [S], to be the vector space of all linear combinations c of differential functionals

∂ k 1 +•••+kn ∂x k 1 1 ...∂x kn n with k 1 + • • • + k n k such that: (a) D 0 q [S] = span({ ∂ 0 ∂x 0 1 ...∂x 0 n }), (b)
c applied to f i , evaluated at q is zero for all integers 1 i m, and (c) for all i ∈ {1, . . . , n}, the anti-differentiation transformation φ j applied to c is in

D k-1 q [S]. Definition 2.2.10 ([DLZ2011, Definition 1]). Let F ∈ C ∞ (R n , R k ) such that F -1 (0)
is a finite set and let a ∈ R n be a solution of the system S = {F = 0}.

Consider the ascending chain of dual spaces

D 0 q [F ] ⊆ D 1 q [F ] ⊆ . . . D h q [F ] ⊆ . . . . If there exists an integer α such that D α q [F ] = D α+1 q [F ],
then the dimension of the vector space D α a [F ] is called the multiplicity of q in the system S. If such an α does not exist, the multiplicity is, by convention, infinity.

As mentioned before, for polynomial systems, the previous definition is equivalent to Definition 2.1.1 [DLZ2011, Theorem 2]. In addition, the following proposition shows that algebraic tools can be used in the smooth case.

Proposition 2.2.11 ([DLZ2011, Corollary 3]). For an integer

k n, let F = (f 1 , . . . , f k ) ∈ C ∞ (R n , R k
) and let a ∈ R n be a solution of the system {F = 0}. Suppose that the multiplicity of a in {F = 0} is m < ∞, then the intersection multiplicity at a of the polynomial system {G = (g 1 , . . . , g k ) = 0} is also m, where g i is equal to the Taylor expansion of f i at a up to degree at least m.

Transversality theorems

The goal of this section is to introduce the concept of genericity and some related results. This has of course a direct relation with our assumptions.

Let U ⊆ R n be an open subset. The set C ∞ (U, R) can be equipped with the so called weak topology that is defined in [Dem2000, §3.9.2] as follows: Let f i be a sequence of elements in C ∞ (U, R). The latter sequence is said to be convergent to a function f if over every compact subset of U , the sequence

∂ k f i ∂x k 1 1 ...x kn n uniformly converges to ∂ k f ∂x k 1 1 ...x kn n , with k = n i=1 k i for all k 1 , . . . , k n in N.
C ∞ (U, R). Naturally, for a positive integer m, the set C ∞ (U, R m ) is equipped with the product topology as C ∞ (U, R m ) can be seen as is the Cartesian product of m copies of C ∞ (U, R).
The key to prove the genericity of our assumptions is Thom's Transversality Theorem. We thus first recall, in this Section, the basics of transversality theory using the notation of Demazure's book [Dem2000].

We work with the set of smooth functions C ∞ (R n , R n-1 ) with the weak (or compact-open) topology [Dem2000, §3.9.2], that is, the convergence is understood as uniform on compact subsets and for any derivative.

A subset of C ∞ (R n , R n-1 ) is called residual if it contains the intersection of a countable family of dense open subsets. The space C ∞ (R n , R n-1 ) is a Baire space [Dem2000, Proposition 3.9.3], that is, every residual subset of C ∞ (R n , R n-1 ) is dense. Definition 2.2.12. A property is generic in C ∞ (R n , R n-1 ) if it is satisfied by a residual subset.
Definition 2.2.13 ([Dem2000, §3.8.3]). Let E R n and F be two finitedimensional real vector spaces and let r 0 be an integer. Let P r (E, F ) be the vector space of polynomial functions of degree at most r from E to F . For an open subset U of E (with respect to the usual topology on E), let J r (U, F ) = U × P r (E, F ) be the space of jets of order r of functions from U to F . Notice that J r (U, F ) can be identified with an open subset of R N for some positive integer N . Let f : U → F be a smooth function, the jet of order r of f is the function

j r f : U ⊂ R n → J r (U, F ) ⊆ R N x → x, f (x), ∂f ∂x 1 (x), . . . , ∂f ∂x n (x), ∂ 2 f ∂x 1 ∂x 2 (x), . . . , ∂ r f ∂x r n (x) .
Let W be a sub-manifold of J r (U, F ). We say that j r f is transverse to W if for all a ∈ U either j r f (a) ∈ W or every vector of R N can be written as a sum of a vector of T j r f (a) W and a vector in the image of the function T a j r f , where T j r f (a) W is the tangent space of W at j r f (a) and T a j r f is the derivative function of j r f at a.

Theorem 2.2.14 (Thom's Transversality Theorem [Dem2000, Theorem 3.9.4]).

Let E and F be two finite-dimensional vector spaces with U an open set in E.

Let r 0 be an integer and W be a sub-manifold of J r (U, F ). Then, the set of

functions f ∈ C ∞ (U, F ) such that j r f is transverse to W is a dense residual subset of C ∞ (U, F ). Proposition 2.2.15 ([Dem2000, Corollary 3.7.3]). Let U be an open subset of R n , N
1 be an integer and W be a sub-manifold of the vector space R N of pure co-dimension m. Assume that the smooth function g :

U → R N is transverse to W , then g -1 (W ) is a (possibly empty) sub-manifold of dimension n -m.
The idea of our proofs of genericity of our assumptions is to express each of them as a system of equations in the jet space. When this system defines a manifold, Thom's theorem applies directly to pull back the manifold from the jet space to the ambient space of the curve and obtain the subset where the assumption is satisfied together with its dimension according to Proposition 2.2.15. A difficulty occurs when the system does not define a manifold. The following corollary overcomes this difficulty in the special case where the system is defined by analytic functions, in other words the system defines an analytic variety. Such a variety does not need to be a manifold but, using the Whitney stratification theorem [Whi1965], the variety is written as a union of manifolds on which Thom's theorem is then applied.

Corollary 2.2.16. Let E and F be two finite-dimensional vector spaces with E of dimension n and U an open set in E. Let r 0 be an integer and W be an analytic variety of J r (U, F ) with co-dimension larger than n, then for a generic P ∈ C ∞ (U, F ), the pre-image of W under j r P is empty.

Proof. Let W = m i=1
W i be a Whitney stratification of W , where the W i 's are submanifolds. Since codim(W ) > n, we have that codim(W i ) > n for any integer

1 i m. Let Γ i = {P ∈ C ∞ (U, F ) | j r P is transverse to W i } and Γ = m i=1 Γ i .
By Theorem 2.2.14, Γ i is residual and so is Γ. Moreover, by Proposition 2.2.15, for P ∈ Γ i the pre-image of W i under j r P is empty. Hence,

(j r P ) -1 (W ) = m i=1 (j r P ) -1 (W i ) = ∅. Lemma 2.2.17. [GG1973, Lemma II.4.3] Let X, Y be smooth manifolds, W ⊂ Y a sub-manifold, and f : X → Y be a smooth function. Let q ∈ X and f (q) ∈ W . Suppose there is an open subset U of f (q) in Y and a submersion φ : U → R k (k = codim(W )) such that W ∩ U = φ -1 (0). Then, f is transversal to W at q iff φ • f is submersion at q.
We will also need a refined version of Thom's theorem in a multijet setting, that is for several points in the source space simultaneously. We give the formal def-initions of the multijet space and function but we do not restate Theorem 2.2.14, Proposition 2.2.15 and Corollary 2.2.16 that also hold for multijets. Definition 2.2.18 ([Dem2000, §3.9.6]). Let U be an open subset of R n and k 1 be an integer. We denote ∆ (k) (U ) the subset of U k consisting of sequences (a 1 , . . . , a k ) of pairwise distinct points of U . For an integer r 0 and a finite dimensional space F , the k-multijet space of order r, J r

(k) (U, F ), is the subset of J r (U, F ) k = (U × P r (E, F )) k consisting of the k-tuples ((a 1 , p 1 ), . . . , (a k , p k )), with (a 1 , . . . , a k ) ∈ ∆ (k) (U ). Let f : U → F be a smooth function, the k-multijet of order r of f is the function j r (k) f : ∆ (k) (U ) → J r (k) (U, F ) (a 1 , . . . , a k ) → (j r f (a 1 ), . . . , j r f (a k )).
Theorem 2.2.19 (Multijet Transversality Theorem [Dem2000, Theorem 3.9.7]). Let E and F be two finite dimensional vector spaces and let U be an open subset of E. Let r 0 and k 1 be two integers and W be a sub-manifold of J r The goal of this thesis is to design certified numerical methods to isolate the singularities of plane curves. Numerical algorithms usually cannot handle singularities. However, under some assumptions, we succeed to handle in a certified way some, although not all, singularities of projections of smooth curves defined in R n .

(k) (U, F ). The set of f ∈ C ∞ (U, F ) such that j r (k) f is transverse to W is a dense residual subset of C ∞ (U, F ).
In this chapter, we introduce a set of assumptions on a curve C in R n , with n > 2 (Section 3.1.2). We study, under these assumptions, the singularities of the plane projection of C classifying them according to Arnold's list [AGZV2012]. This classification is later used to show a correspondence between the singularities of the plane projection of C and the solutions of a regular and square system in R 2n-1 introduced in Chapter 4. In section 3.2.1, we prove that the assumptions mentioned above are satisfied by a generic curve in R n . In addition, in Section 3.2.2, we prove that the silhouette of a generic surface in R n satisfies a part of these assumptions. we conclude the chapter by a conjecture that the silhouette of a generic surface satisfies the rest of the assumptions.

Notation and assumptions

in R n . Let C ∞ (R n , R n-1
) denote the set of smooth functions (i.e., differentiable infinitely many times) from R n to R n-1 . Consider the function

P = (P 1 , . . . , P n-1 ) ∈ C ∞ (R n , R n-1
). We will denote by C (resp. C) the solution set of the system {P 1

(x) = • • • = P n-1 (x) = 0}, with x = (x 1 , . . . , x n ) ∈ B (resp. with x ∈ B). Also, consider the projection π C (resp. π C ) from C (resp. C) to the (x 1 , x 2 )-plane. Unless otherwise stated, the plane projection of a point x ∈ R n is (x 1 , x 2 ).
If C is a smooth curve (see the definition below), define L c (resp. L c ) to be the set of points q in C (resp. C) such that the tangent line, denoted by T q C, (resp. T q C) is orthogonal to the (x 1 , x 2 )-plane. We also define the set L n (resp. L n ) to be the set of points q in C (resp. C) such that the cardinality of the pre-image of π C (q) under π C (resp. π C ) is at least two. We will see later that, under some generic assumption, L c (resp. L n ) is equal to the set of points in C that project to a cusp (resp. node) which justifies the subscript c (resp. n).

Assumptions

Recall that our goal is to design a numerical algorithm for isolating the singularities that appear in the plane projection of a curve C in R n . Numerical algorithms usually cannot handle degenerate cases, that is, singularities in our context. However, under some assumptions on C, we succeed to isolate in a certified way some singularities of the projection. Namely, we require that the singularities are "generic", that is, only nodes can appear in the projection (Assumption A 5 ). Our other assumptions on C are, roughly speaking, that it is smooth (A 1 ), that its projection only has a discrete set of singularities (A 4 ), and that at most two points of C project on each singularity (A 3 ). We will prove in Chapter 5 that our numerical algorithm is certified and terminates under these assumptions and in Section 3.2 that these assumptions are generically satisfied.

Recall that we denote by J P the Jacobian matrix of the function P in C ∞ (R n , R n-1 ). Consider the following assumptions:

A 1 -For all q ∈ C, rank(J P (q)) = n -1. In particular, C is a smooth curve. 5A 2 -The set L c is discrete and does not intersect the boundary of B.

A 3 -For all points p = (α, β) ∈ π C (C), the pre-image of p under π C consists of at most two points in B counted with multiplicities in the system

{P (x) = 0 ∈ R n-1 , x 1 -α = x 2 -β = 0}.
A 4 -The set L n is discrete and does not intersect the boundary of B.

A 5 -The singular points of π C (C) are only nodes (see Definition 2.2.4).

We also consider a weaker version of Assumption A 5 in which ordinary cusps can also appear in the projection:

A - 5 -
The singular points of π C (C) are only ordinary cusps or nodes (see Definition 2.2.4).

Definition 3.1.1. Assumptions A 1-5 are called the strong assumptions and Assumptions A 1-4 and A - 5 are called the weak assumptions.

Figure 3.1 illustrates the assumptions that we have just introduced. Our goal in the next section is to prove that a generic curve C satisfies our assumptions. Loosely speaking, the genericity implies that a "random" curve satisfies these assumptions. Moreover, if a curve does not satisfy the assumptions, there exists a sequence of generic curves that is converging to that curve.

Remark 3.1.2. Our motivation for also considering these weak assumptions is dual. First, our certified algorithm for isolating the singularities of the projection of curves satisfying the strong assumptions also works, to some extent, if only the weak assumptions hold: namely, it outputs a superset of the isolating boxes of the singularities. Second, we conjecture that our weak assumptions are satisfied by silhouette curves of generic surfaces (see Proposition 3.2.17 and Conjecture 3.2.18).

Genericity of the assumptions

Recall from the introduction (namely in Section 1.2.2) that a property in

C ∞ (R n , R n-2 ) is generic if it is satisfied by all elements of a residual subset (a countable intersection of open dense subsets) in C ∞ (R n , R n-2 ).
The key to prove the genericity of our assumptions mentioned on Section 3.1.2 is Thom's Transversality Theorem 2.2.14. We mentioned, in Section 2.2.4, the basics of transversality theory using the notation of Demazure's book [Dem2000]. We prove, in Section 3.2.1, that all assumptions of Section 3.1.2 are satisfied for a generic curve. In Section 3.2.2, we consider the special case where the curve is the silhouette of

q points in L n π C (C) point in L c R n-2 x1 x2
(a) Assumption A1: C is a smooth; the rank of the Jacobian of P at q is full. a surface and prove that Assumptions A 1 , A 2 , A 4 are generically satisfied in this case; the genericity of Assumptions A 3 and A - 5 is left as a conjecture.

Genericity of the assumptions for a curve in R n

We are going to prove that each assumption in Section 3.1 is generic. Hence, the combination of these assumptions is also generic since a countable intersection of residual subsets in

C ∞ (R n , R n-1 ) is residual. Lemma 3.2.1. Assumption A 1 is generic.
Proof. Consider the jet of order 1 of the function P ∈ C ∞ (R n , R n-1 ):

j 1 P : R n → J 1 (R n , R n-1 ) = R n × R n-1 × R (n-1)×n
x → (x, P (x), J P (x)) = (x, y, z).

We represent the jet space by the variables x ∈ R n , y ∈ R n-1 and z ∈ R (n-1)×n . With abuse of notation, we can see the variable z as a (n -1) × n-matrix. Define the variety

W = {(x, y, z) ∈ R n × R n-1 × R (n-1)×n | y = 0, rank(z) n -2}.
The variety W is a product of a determinantal variety in R (n-1)×n of dimension n 2 -n -2 (by Proposition 2.1.4) and a linear space of dimension n in R n × R n-1 . Thus, W is a variety of co-dimension n + 1 in R n × R n-1 × R (n-1)×n . Hence, by Corollary 2.2.16, there exists a residual subset Γ 1 ⊂ C ∞ (R n , R n-1 ), such that for P ∈ Γ 1 the pre-image of W under j 1 P is empty. Consequently, for a generic P ∈ C ∞ (R n , R n-1 ) and any q ∈ C, we have that q / ∈ (j 1 P ) -1 (W ) = ∅, thus rank(J P (q)) = n -1, which is Assumption A 1 .

The following lemma proves a stronger claim than the genericity of Assumption A 2 , namely, that generically, L c is empty. However, as mentioned before, being L c discrete, is the case for the Silhouette of a generic surface in R n (more precisely, generic P in C ∞ (R n , R n-2 )) which justify why we choose the weaker assumption. Lemma 3.2.2. Assumption A 2 is generic. Moreover, generically, the set L c is empty.

Proof. We consider the jet of order 1 of the function P ∈ C ∞ (R n , R n-1 ) as is the proof of Lemma 3.2.1 with the same notation. Define the matrix T 1 (z) (resp. T 2 (z)) to be the sub-matrix of z obtained by removing the first (resp. second) column. Consider the variety W ⊂ J 1 (R n , R n-1 ) defined by {y = 0 ∈ R n-1 , det(T 1 (z)) = det(T 2 (z)) = 0}. Notice that L c is included in the pre-image of W under j 1 P since L c is the set of points of the curve C that are both x 1 and x 2 -critical. By Lemma 2.1.5, we have that both det(T 1 (z)) and det(T 2 (z)) are irreducible polynomials. By [CLO2007, §9.4 Prop 10], a proper sub-variety of an irreducible variety is of lower dimension, we deduce that the common zero locus of det(T 1 (z)) and det(T 2 (z)) is of co-dimension at least two. We deduce that codim(W ) > n. By Corollary 2.2.16, there exists a residual subset Γ 2 ⊂ C ∞ (R n , R n-1 ), such that for P ∈ Γ 2 ∩Γ 1 , the pre-image of W under j 1 P is empty and hence L c is empty, which implies Assumption A 2 .

The following lemma provides a characterization of the points in L c (under Assumption A 1 ) in terms of multiplicity. Such an equivalence is useful to ease proofs and reduce technicalities.

Lemma 3.2.3. Let P = (P 1 . . . , P n-1 ) ∈ C ∞ (R n , R n-1 ) satisfy Assumption A 1 . Let q be in C such that the multiplicity of the system S = {P (x) = 0 ∈ R n-1 , x 1 - α = x 2 -β = 0} at q is finite, where (α, β) = π C (q) ∈ R 2 .
Then, q ∈ L c if and only if the multiplicity of the system S at q is at least two.

Proof. Without loss of generality assume that q = 0 ∈ R n .

Sufficiency: Assume that q ∈ L c . Let v = (v 1 , . . . , v n ) be a non-trivial vector of the tangent line of C at q. Thus, J P (q) • v T = 0. By the definition of L c we have v 1 = v 2 = 0. Define the differential operator c = n i=3 v i ∂ ∂x i . Notice that c • P j = n i=3 v i ∂P j ∂x i (q) = 0 for all integers 1 j n -1 (see [DLZ2011, 2.1] for
the definition of c • P j ). Moreover, by the definition of c and since

v 1 = v 2 = 0, we have c • (x 1 ) = c • (x 2 ) = 0. Hence, c ∈ D 1 q [S] \ D 0 q [S]
. Thus, dim(D 1 q ) > 1. Hence, the multiplicity of S at q is at least two.

Necessity: Assume that the multiplicity of S at q is at least two, then

D 0 q [S] D 1 q [S]
. This implies that there exists a non-trivial differential operator c

= n i=1 c i ∂ ∂x i ∈ D 1 q [S] \ D 0 q [S] such that:
(a) We have that c • P j = 0 for all integers 1 j n -1 which implies that if we write v i = c i , with 1 i n, the non-trivial vector v is in the tangent space of C at q.

(b) We have that c • (x 1 ) = c • (x 2 ) = 0, equivalently, c 1 = c 2 = 0. Thus, v 1 = v 2 = 0.
The tangent line to the curve at q is thus orthogonal to the (x 1 , x 2 )-plane. Thus, q ∈ L c . Lemma 3.2.4. Assumption A 3 is generic.

Proof. Let us consider the 3-multijet of order 0:

j 0 (3) P : ∆ (3) (R n ) → J 0 (3) (R n , R n-1 ) = (R n × R n-1
) 3 (x, x , x ) → ((x, P (x)), (x , P (x )), (x , P (x ))) = ((x, y), (x , y ), (x , y ))

where every element in the jet space J 0

(3) (R n , R n-1 ) is of the form ((x, y), (x , y ), (x , y )), where x = (x 1 , . . . , x n ), x , x ∈ R n and y, y , y ∈ R n-1 . Consider the linear sub-manifold W = {x 1 = x 1 = x 1 , x 2 = x 2 = x 2 , y = y = y = 0}, the co-dimension of W is thus 3n+1 which is larger than the dimension of the source space ∆ (3) (R n ) which is 3n. Thus, by Corollary 2.2.16, there exists a residual subset Γ 3 ⊂ C ∞ (R n , R n-1 ), such that for P ∈ Γ 3 , the pre-image of W by j 0

(3) is empty, which translates to the fact that there are no pairwise distinct points q, q , q in C such that π C (q) = π C (q ) = π C (q ). This is also equivalent to say that the system

S = {P (x) = 0 ∈ R n-1 , x 1 -α = x 2 -β = 0} has at most two distinct solutions (without counting multiplicities) for any (α, β) ∈ R 2 .
Using Γ 1 , Γ 2 as defined in the proofs of Lemmas 3.2.1 & 3.2.2 and Γ 3 defined above, we define Γ 4 = Γ 1 ∩ Γ 2 ∩ Γ 3 which is thus a residual set and let P be in Γ 4 . Since P is in Γ 3 , the system S has at most two distinct solutions. In addition, since P is in Γ 2 ∩ Γ 1 , one has that L c is empty and finally together with Lemma 3.2.3, since P is in Γ 1 , this implies that these solutions have multiplicity exactly 1 is S. For P in the residual set Γ 4 , the number of solutions counted with multiplicities of S is thus at most 2, which is Assumption A 3 . Lemma 3.2.5. Assumption A 4 is generic.

Proof. Let us consider the 2-multijet of order 0 of P :

j 0 (2) P : ∆ (2) (R n ) → J 0 (2) (R n , R n-1 ) = (R n × R n-1 ) 2 (x, x ) → ((x, P (x)), (x , P (x ))) = ((x, y), (x , y ))
where every element in the jet space J 0

(2) (R n , R n-1 ) is of the form ((x, y), (x , y )),

where x = (x 1 , . . . , x n ), x ∈ R n and y, y ∈ R n-1 . Consider the linear sub- manifold W = {x 1 = x 1 , x 2 = x 2 , y = y = 0} of the jet space J 0 (2) (R n , R n-1 ). Notice that, (j 0 (2) P ) -1 (W ) contains the set L n = {(q 1 , q 2 ) ∈ ∆ (2) (R n ) ∩ C × C | π C (q 1 ) = π C (q 2 )
} and L n is the image of L n by the projection (q 1 , q 2 ) → q 1 . We have dim(∆ (2) (R n )) = 2n and, since W is linear, its co-dimension is easily computed codim(W ) = 2(2n -1) -(2 + 2(n -1)) = 2n. Proposition 2.2.15 thus yields that generically (j 0

(2) P ) -1 (W ) is a sub-manifold of dimension zero that is a discrete set in R n , and so is L n . Now, we prove that, generically, L n does not intersect the boundary of B. The boundary ∂B of the box B is included in the union of the supporting hyperplanes

H i of its 2 n faces of dimension n -1, that is ∂B = ∪ i=2 n i=1 H i . Define the linear sub-manifold W i = {((x, y), (x , y )) ∈ W | x ∈ H i or x ∈ H i },
notice that this adds one equation to W and thus increases the co-dimension of W by one, thus codim(W i ) = 2n + 1. By Corollary 2.2.16, we have that, generically, the preimage of W i under j 0

(2) P is empty, which translates to the fact that there is no point of L n on ∂B ∩ H i . This is also true for any i and thus, generically, L n does not intersect the boundary of B.

Proving the genericity of Assumptions A 5 and A - 5 is a bit more complicated. We first study the singularity types that occur on the plane curve π C (C) under Assumptions A 1 , A 2 , A 3 and A 4 in Lemmas 3.2.6 and 3.2.7. Then, we translate their results into algebra (Lemma 3.2.9) which paves the way to use transversality theorem in proving the genericity of Assumptions A 5 and A - 5 in Corollary 3.2.10. Lemma 3.2.6. Under Assumptions A 1 , A 2 , A 3 and A 4 , let q ∈ C and p = π C (q). If q ∈ L c ∪ L n , then p is a smooth point of the plane curve π C (C).

Proof. As illustrated in Figure 3.2, since q ∈ L c , the plane projection of T q C is a line, or equivalently, the derivative T q π C of π C at q is injective. Thus, π C is an immersion at q. Hence, for a small enough neighbourhood U 0 of q in R n , we have that π C restricted to V = U 0 ∩ C is embedding (Proposition 2.2.7). We are going to prove that, assuming that U 0 is small enough, the curve π C (C) has exactly one branch around π C (q) which implies that π C (C) is smooth at π C (q) since C is smooth at q by Assumption A 1 .

To prove this claim, assume that there exists an open subset U 0 in R n such that the set V = U 0 ∩ C and V are disjoint, but π C (q) is in the closure of π C (V ). Let q k be a sequence of points in V such that π C (q k ) converges to π C (q). Since B is compact, there exists a convergent sub-sequence of q k that has a limit q in B. Notice that π C (q ) = π C (q) by the continuity of π C . Hence, q, q are both in L n . However, since q ∈ L n , we must have that q ∈ B. Hence, q is in the boundary of B which contradicts Assumption A 4 . Hence, the curve π C (C) has exactly one smooth branch around π C (q) which concludes the proof.

Regarding our goal about the singularities of π C (C), by Lemma 3.2.6, we can concentrate our attention to the plane projections of the points in

L c ∪ L n . Lemma 3.2.7. Under Assumptions A 1 , A 2 , A 3 and A 4 , if q ∈ L n , then π C (q)
is a singular point of the plane curve π C (C). More precisely, either π C (q) is of type A - 2k+1 with k 0, or there exists a non-null smooth function g defined in a neighbourhood of 0 ∈ R with ord(g) = ∞ such that (π C (C), π C (q)) is equivalent to the curve defined by x 2 -g(y 2 ) = 0 at the origin. Proof. As in Figure 3.4, let p = π C (q), according to A 3 , π -1 C (p) has at most two points, and since q is in L n , it also has at least two points. Define q such that π -1 C (p) = {q, q } and denote the plane curve π C (C) by C. Without loss of generality, one can assume that p = (0, 0). In addition, A 3 also implies that the multiplicities of q and q in the system {P (x) = 0 ∈ R n-1 , x 1 = x 2 = 0} are one. With Assumption A 1 , Lemma 3.2.3 then implies that the tangents to C at q and q are not orthogonal to the (x 1 , x 2 )-plane. Thus there exists two neighborhoods N q and N q of q and q in R n such that π restricted to C∩N q (resp. C∩N q ) is an embedding. Let D k be a sequence of open disks centered at p and of radius 1 k . By contradiction, if for all k, there exists points q k ∈ C such that q k is not in N q ∪ N q and π C (q k ) ∈ D k , then the limit q ∞ is a point of C distinct from q and q , and π C (q ∞ ) = p. If q ∞ is in B, it contradicts A 3 and if it is in B, it contradicts A 4 . Thus for a small enough neighborhood of p, the projection of the curve is restricted to the projection of the two branches around q and q . Finally, if for all D k , the pre-image of π -1 (D k ) contains a point in L n \ {q, q }, then this contradicts the discreteness assumption A 4 . Thus there exists a neighborhood N ⊆ R 2 of p such that π -1 C (N ) is a union of two smooth (Assumption A 1 ) open subsets of C such that q is on one branch and q on the other, and π C restricted to π -1 C (N ) \ {p, q } is an embedding. The projection of these two smooth branches are thus two smooth curves in the plane. Let these two smooth plane branches be defined by the zero sets of the smooth functions f 1 and f 2 in C ∞ (R 2 , R). Let u (resp. u ) be a non-zero tangent vector of C at q (resp. q ) and v (resp. v ) be its projection in R 2 . We distinguish two cases: (a) The vectors v and v are independent in R 2 (see Figure 3.3). Thus, v and v give rise to a local coordinate system (x, y) in a neighborhood of p in R 2 . The vector v being tangent to the zero set of f 1 , one has ∂f 1 ∂x (p) = 0 and ∂f 1 ∂y (p) = 0. By the implicit function theorem [Dem2000, Corollary 2.7.3.], we deduce that there exists a real smooth function h 1 such that y = x 2 • h 1 (x) is a local parameterization of the zero set of f 1 . Similarly, there exists a smooth function h 2 such that x = y 2 • h 2 (y) is a local parameterization of the zero set of f 2 . Thus (x, y)

U 0 T q C q R n V B R 2 p T p π C (C)
∈ N iff f (x, y) = f 1 (x, y)f 2 (x, y) = 0 iff (y -x 2 • h 1 (x))(x -y 2 • h 2 (y)) = 0, equivalently, [y -x -x 2 • h 1 (x) + y 2 • h 2 (y)] 2 -[y + x -x 2 • h 1 (x) -y 2 • h 2 (y)] 2 = 0. The change of coordinates X = y -x + x 2 • h 1 (x) + y 2 • h 2 (y) and Y = y + x + x 2 • h 1 (x) -y 2 • h 2 (y) is a diffeomorphism since det(J x,y (X, Y )) p = 0.
Then, the local equation of the curve C at p is of the form X 2 -Y 2 with these new coordinates, which means that p is a A - 1 or node singularity.

(b) The vectors v and v are co-linear (see Figure 3.4). Then, choose v ∈ T p R 2 linearly independent from v, the vectors v, v give rise to a coordinate system (x, y) at p. In this coordinate system, we thus have

∂f 1 ∂x (p) = ∂f 2 ∂x (p) = 0, ∂f 1
∂y (p) = 0 and ∂f 2 ∂y (p) = 0. By the implicit function theorem, there exist smooth functions h 1 and h 2 such that locally f (x, y) = 0 if and only if (y -

x 2 • h 1 (x))(y -x 2 • h 2 (x)) = 0. The last equality is equivalent to (2y - x 2 (h 1 (x) + h 2 (x))) 2 -x 4 (h 1 (x) -h 2 (x)) 2 = 0.
Assumption A 4 ensures that the projections of the two branches have only one common point, such that h 1 (x) -h 2 (x) does not vanish identically. We distinguish two cases:

q p N u q u v = k • v v
(i) mult(h 1 (x) -h 2 (x)) = k ∞, then h 1 (x) -h 2 (x) = x k • u with u(p) = 0 and without loss of generality, assume that u(p) > 0. The change of coordinates X = 2y -x 2 (h 1 (x) + h 2 (x)) and Y = x • u 1 2+k is a diffeomorphism (notice that indeed u 1 2+k is a smooth function around p).
Then, the local equation of the curve C at p is of the form X 2 -Y (2k+3)+1 with these new coordinates, which means that p is a singularity of type

A - 2k+3 . (ii) mult(h 1 (x)-h 2 (x)) = ∞. Since the function x 4 (h 1 (x)-h 2 (x)) 2 is even,
by Theorem 2.2.2, there exists a smooth function g such that x 4 (h 1 (x)h 2 (x)) 2 = g(x 2 ). Thus, taking the diffeomorphism X = 2y -x 2 (h 1 (x) + h 2 (x)) and Y = x, we get the second case of the claim.

The next definition and lemma are technical tools for proving the genericity of Assumptions A - 5 and A 5 . We want to translate the result of Lemma 3.2.7 to the language of algebra. The following definition is needed for this characterization: Definition 3.2.8. Consider P = (P 1 , . . . , P n-1 ) ∈ C ∞ (R n , R n-1 ) satisfying Assumption A 1 and recall that we denote by J P (q) the Jacobian matrix of P at the point q. We define the (n -1) × (n -2) sub-matrix M 1,2 (q) obtained by removing the first two columns of J P (q) and the (n -1) × 2 sub-matrix N P (q) formed by the first two columns of J P (q). Let q 1 , q 2 ∈ C, we define the square matrix of size 2n -2, M (q 1 , q 2 ) = N P (q 1 ) 0 M 1,2 (q 1 ) N P (q 2 ) M 1,2 (q 2 ) 0 .

Lemma 3.2.9. Using the same assumption and notation as in Definition 3.2.8, let q 1 and q 2 be distinct points of C with π C (q 1 ) = π C (q 2 ), then M (q 1 , q 2 ) is invertible if and only if none of q 1 or q 2 is in L c and the plane projections of the tangent lines of C at q 1 and q 2 do not coincide.

Proof. We prove the converse statement using det(M (q 1 , q 2 )) = 0 ⇐⇒ There exist α ∈ R 2 and β, γ ∈ R n-2 such that the vector x = (α, β, γ) is not trivial and satisfies the equality M (q 1 , q 2 ) • x T = 0. ⇐⇒ (α, β) and (α, γ) are in the tangent lines T q 1 C and T q 2 C respectively and at least one of them is not trivial. The last statement can be split in two cases:

• α is not trivial which is equivalent to say that the plane projections of T q 1 C and T q 2 C are both generated by α and coincide.

• α = (0, 0) which is equivalent to β or γ is not trivial, which is equivalent to T q 2 C or T q 1 C projects to a point in the plane, which is equivalent to q 1 or q 2 is in L c . Now, we are ready to prove the genericity of Assumptions A 5 -and A - 5 -:

Corollary 3.2.10. Assumption A 5 is generic. Consequently, so is Assumption A - 5 .

Proof. Let B be an open n-box. Recall that generically L c (and hence L c ) is empty (Lemma 3.2.2). Hence, it is enough to prove that for a generic P ∈ C ∞ (R n , R n-1 ), the singular points of π C (C) are only nodes (recall that by Lemma 3.2.6, under the generic assumptions

A 1 , A 2 ,A 3 and A 4 , the points in C \ (L c ∪ L n ) project to smooth points). Let Γ 0 be the set of P ∈ C ∞ (R n , R n-1
) such that P satisfies Assumptions A 1 , A 2 , A 3 and A 4 . The previous lemmas of this section show that Γ 0 is residual in

C ∞ (R n , R n-1 ).
Let us consider the 2-multijet of order 1 of P :

j 1 (2) P : ∆ (2) (R n ) → J 1 (2) (R n , R n-1 ) ⊆ (R n × R n-1 × R (n-1
)×n ) 2 (x, x ) → ((x, P (x), J P (x)), (x , P (x ), J P (x ))) = ((x, y, z), (x , y , z )) Let s, s (resp. r, r ) be the sub-matrices of z, z respectively obtained by removing the first two columns (resp. obtained by the first two columns). Define the matrix M = r 0 s r s 0 and the variety

W = {((x, y, z), (x , y , z )) ∈ (R n ×R n-1 ×R (n-1)×n ) 2 | y = y = 0, x1 = x 1 , x2 = x 2 , det(M ) = 0}.
The variety W is a product of a determinantal variety and a linear space, thus its co-dimension is

codim(W ) 2n + 1 > 2n = dim(∆ (2) (R n ))
. Hence, by Corollary 2.2.16, there exists a residual subset Γ 0 in C ∞ (R n , R n-1 ) such that for all P ∈ Γ 0 , the pre-image of W under j 1

(2) P is empty. Let then P be in the residual set Γ 0 ∩ Γ 0 . By Lemma 3.2.9 and since L c is empty, we deduce that for distinct q 1 , q 2 ∈ C with π C (q 1 ) = π C (q 2 ), the plane projections of the lines T q 1 C and T q 2 C intersect transversely if and only if j 1

(2) ((q 1 , q 2 )) ∈ W . Finally, by Lemma 3.2.7 (Step (a) of the proof), we deduce that π C (q 1 ) = π C (q 2 ) is a node in π C (C).

Genericity of the assumptions for the silhouette of a surface in R n

In this section, we focus on the special case of silhouette curves of surfaces in R n . For an open n-box B and P in C ∞ (R n , R n-2 ) such that S = P -1 (0) is a smooth 2-sub-manifold in R n , the silhouette of P is the set of points q of this surface S such that the projection (with respect to a fixed direction) of the tangent plane T q S to R 2 is not surjective (see Figure 3.5). We prove that Assumptions A 1 , A 2 & A 4 are satisfied for a generic silhouette, and we only conjecture that Assumptions A 3 & A - 5 also hold generically. We start by formalizing algebraically the definition of the silhouette curve. Definition 3.2.11. For an integer n 3, let P = (P 1 , . . .

, P n-2 ) ∈ C ∞ (R n , R n-2
). Define the smooth function P n-1 = det ∂P i ∂x j 1 i n-2 3 j n and P = (P 1 , . . . , P n-1 ). We define the curve C (and C) as in Section 3.1 and call it the silhouette of P .

As mentioned before, we want now to prove that the silhouette of a generic surface in R n , Assumptions A 1 , A 2 and A 4 are satisfied. The argument used for the proofs are more involved and complicated comparing with the ones presented in Section 3.2.1. For Assumption A 1 , the basic idea of the proof of Lemma 3.2.12 is Lemma 2.2.17. Proposition 3.2.12. For a generic P ∈ C ∞ (R n , R n-2 ), the function P satisfies Assumption A 1 . Proof. Consider the jet of order 1 of P :

j 1 P : R n → J 1 (R n , R n-2 ) = R n × R n-2 × R (n-2)×n R n 2 -2 = R N
x → (x, P (x), J P (x)) = (x, y, z).

We represent the jet space by the vectors x ∈ R n , y ∈ R n-2 and the ((n -2) × n)matrix z ∈ R (n-2)×n . Let T (z) denote the sub-matrix obtained by removing the first two columns of z. Define the variety W = {y = 0, det(T (z)) = 0} = {y = 0, rank(T (z))

n -3} in R N . According to Proposition 2.1.4, W = Reg(W ) ∪ Sing(W ) where Reg(W ) (resp. Sing(W )) is the set of smooth (resp. singular) points in W and

Reg(W ) = {(x, y, z) ∈ R N | y = 0, rank(T (z)) = n -3} Sing(W ) = {(x, y, z) ∈ R N | y = 0, rank(T (z)) < n -3}.
In addition, Proposition 2.1.4 yields that Reg(W ) is a manifold of co-dimension n -1 and Sing(W ) is a variety of co-dimension n + 2. Since the co-dimension of Sing(W ) is larger than that of the source space, Corollary 2.2.16 implies that, generically, (j 1 P ) -1 (Sing(W )) = ∅. One thus has (j 1 P ) -1 (W ) = (j 1 P ) -1 (Reg(W )). Consider the function

ϕ : R n × R n-2 × R (n-2)×n → R n-2 × R χ = (x, y, z) → (y, det(T (z))), such that ϕ -1 (0) = W . Its Jacobian matrix is J ϕ = 0 (n-2)×n I (n-2)×(n-2) 0 (n-2)×(n-2)n 0 1×(n) 0 1×(n-2) v(z)
, where 0 k 1 ×k 2 (resp.

I k 1 ×k 2 ) is
the zero (resp. identity) matrix of size k 1 × k 2 and the vector v(z) is the adjugate matrix of T (z) written as the concatenation of its lines:

v(z) = (Adj ij (T (z))) 1 i n-2 3 j n ∈ R (n-2) 2 . Let χ = (x, y, z) ∈ Reg(W ),
then rank(T (z)) = n -3, thus there exists a pair (i, j) such that Adj ij (T (z)) = 0. Hence, the vector v(z) is non-trivial and J ϕ (χ) has full rank n -1. The function ϕ is thus a submersion (Definition 2.2.6) on Reg(W ). Theorem 2.2.14 yields that, generically, j 1 P is transverse to the manifold Reg(W ). As illustrated in Figure 3.6, since ϕ is a submersion on Reg(W ), by Lemma 2.2.17, we imply that that P = ϕ • j 1 P is a submersion on (j

1 P ) -1 (Reg(W )) = (j 1 P ) -1 (W ) = (j 1 P ) -1 (ϕ -1 (0)) = (ϕ • j 1 P ) -1 (0) = P -1 (0) = C.
In other words, J P has full rank n -1 on C, which is Assumption A 1 .

Proposition 3.2.13. For a generic P ∈ C ∞ (R n , R n-2 ), the function

P satisfies Assumption A 2 . R n-2 × R P = ϕ • j 1 P R N = R n × R n-2 × R (n-2)×n R n j 1 P W ϕ Figure 3
.6: Illustration of Lemma 2.2.17 which is used in the proof of Proposition 3.2.12. Since ϕ is a submersion over Reg(W ) and since j 1 P is, generically, transversal to Reg(W ), the lemma implies that the composition P = ϕ • j 1 P is generically a submersion.

Proof. First we prove that, generically, L c is discrete. For any

P ∈ C ∞ (R n , R n-2 ) consider j 2 P : R n → J 2 (R n , R n-2 ) ⊂ R n × R n-2 × R (n-2)×n × R n 2 (n-2) = R N . Assume that every element in R N is represented as (x, y, z, h), where x ∈ R n , y ∈ R n-2 , z ∈ R (n-2)×n and h ∈ R n 2 (n-2) .
With abuse of notation we can consider z as a ((n -2) × n)-matrix. Let T (z) denote the matrix obtained by removing the first two columns of z. The Jacobian matrix J P is a function of the derivatives

( ∂P i ∂x j , ∂ 2 P l ∂x k ∂xs ) 1 i,l n-2 1 j,k,s n
, it can thus be seen in the jet space as a function of z and h, J P (z, h). Define the matrix T 1 (z, h) (resp. T 2 (z, h)) to be the submatrix of J P (z, h) obtained by removing the first (resp. second) column. Define the variety

W = {(x, y, z, h) | y = 0 ∈ R n-2 , det(T (z)) = det(T 1 (z, h)) = det(T 2 (z, h)) = 0}, so that L c is included in the pre-image of W under j 2 P . Let W 1 = {(x, y, z, h) | y = 0 ∈ R n-2
, det(T (z)) = 0}, we already showed in the proof of Proposition 3.2.12 that W 1 is an irreducible variety of co-dimension n -1.

In addition, det(T 1 (z, h)) does not identically vanish on W 1 , thus W is a proper sub-variety of the irreducible variety W 1 and [CLO2007, §9.4 Prop 10] implies that codim(W ) > codim(W 1 ) = n -1.

Now, write W = Reg(W ) ∪ Sing(W ),
⊂ C ∞ (R n , R n-2 ) such that if P ∈ Γ , then the pre-image of Sing(W ) under j 2 P is empty. De- fine Γ = { P ∈ C ∞ (R n , R n-2 ) | j 2 P is transverse to Reg(W )} ∩ Γ . Notice that if P ∈ Γ, then L c is included in the pre-image of Reg(W ) under j 2 P . Hence, since codim(Reg(W )) = codim(W )
n, we have by Proposition 2.2.15 that L c is a sub-manifold of dimension, at most, zero. Thus, L c is discrete for all P ∈ Γ. Using Theorem 2.2.14 we deduce that Γ is residual.

The proof that L c does not intersect the boundary of B can be done analogously as in the proof of Lemma 3.2.5. Proposition 3.2.14. For a generic P ∈ C ∞ (R n , R n-2 ), the function P satisfies Assumption A 4 .

Proof. Consider the 2-multijet j 1

(2) P :

∆ (2) (R n ) → J 1 (2) (R n , R n-2 ) = (R n × R n-2 × R (n-2)×n ) 2 of the function P ∈ C ∞ (R n , R n-2 ), where (R n × R (n-2)×n × R (n-2)×n
) 2 is described by the coordinates x, x ∈ R n , y, y ∈ R n-2 and z, z ∈ R (n-2)×n . With abuse of notation we can consider z and z as ((n -2) × n)matrices. Let T (z) (resp. T (z )) denote the matrix obtained by removing the first two columns of z (resp. z ). Define the variety W to be the solution set of the system {y = y = 0, x 1 -

x 1 = x 2 -x 2 = det(T (z)) = det(T (z )) = 0}.
Denote by Reg(W ) the regular part of W . By Proposition 2.1.4 (a) we deduce that W is of co-dimension 2n. Using the same argument in the proof of Proposition 3.2.12, we deduce that there exists a residual set

Γ ⊂ C ∞ (R n , R n-2 ) such that if P ∈ Γ, then the image of ∆ 2 (R n ) under j 1
(2) P is contained in Reg(W ). Moreover, by Proposition 2.2.15, we have that M 1,2 = (j 1

(2) P ) -1 (Reg(W )) = (j 1 (2) P ) -1 (W ) is a sub-manifold of dimension zero in ∆ 2 (R n ). Notice that L n is the image of M 1,2 under the projection (x, x ) → x. Since M 1,2 is of dimension zero, then so is L n . Thus we have just proven that, if P ∈ Γ, then L n is a sub-manifold of dimension zero. Hence, L n is discrete.

The proof that L n does not intersect the boundary of B can be done analogously as in the proof of Lemma 3.2.5.

Assumption A 3 can be rephrased by the three following assumptions:

A 3(a) -There are no pairwise distinct q, q , q ∈ C such that π C (q) = π C (q ) = π C (q ).

A 3(b) -L c ∩ L n = ∅.
A 3(c) -For q ∈ L c , the multiplicity of the system {P (x) = 0 ∈ R n-1 , (x 1 , x 2 ) = π C (q)} at q is exactly two.

Using this rephrasing, we next show that Assumptions A 3(a) & A 3(b) generically hold and we leave the genericity of Assumption A 3(c) as a conjecture.

Proposition 3.2.15. For a generic function

P ∈ C ∞ (R n , R n-2 ), Assump- tion A 3(a) -holds. Proof. Consider the 3-multijet j 1 (3) P : ∆ (3) (R n ) → J 1 (3) (R n , R n-2 ) = (R n × R n-2 × R (n-2)×n ) 3 . Assume that every element in (R n × R n-2 × R (n-2)n
) 3 is of the form ((x, y, z), (x , y , z ), (x , y , z )), where x, x , x ∈ R n , y, y , y ∈ R n-2 and z, z , z ∈ R (n-2)×n . With abuse of notation we can consider z, z and z as ((n -2) × n)-matrices. Let T (z), T (z ), T (z ) denote the matrices obtained by removing the first two columns of z, z , z respectively. Consider the variety W defined by the equations:

{x 1 = x 1 = x 1 , x 2 = x 2 = x 2 , y = y = y = 0 ∈ R n-2 , det(T (z)) = det(T (z )) = det(T (z )) = 0}.
Notice that dim(∆ (3) (R n )) = 3n < 3n + 1 = codim(W ). Hence, by Corollary 2.2.16, we have that, generically, the pre-image of W under j 1

(3) P is empty. Hence, there are no pairwise different q, q , q ∈ C such that π C (q) = π C (q ) = π C (q ). Proposition 3.2.16. For a generic function

P ∈ C ∞ (R n , R n-2 ), Assump- tion A 3(b) -holds. Proof. Consider the 2-multijet j 2 (2) P : ∆ (2) (R n ) → J 2 (2) (R n , R n-2 ) = (R n × R n-2 × R (n-2)×n × R n 2 (n-2) ) 2 of the function P ∈ C ∞ (R n , R n-2 ), where (R n × R (n-2)×n × R (n-2)×n × R n 2 (n-2) ) 2 is described by the coordinates x, x ∈ R n , y, y ∈ R n-2 , z, z ∈ R (n-2)×n and h, h ∈ R n 2 (n-2) .
With abuse of notation we can consider z and z as ((n -2) × n)matrices. Let T (z) (resp. T (z )) denote the matrix obtained by removing the first two columns of z (resp. z ). Define the matrices T 1 (z, h), T 2 (z, h) as in the proof of Lemma 3.2.13 and the variety W to be the solution set of the system {y = y = 0 ∈ R n-2 , x 1 -

x 1 = x 2 -x 2 = det(T (z)) = det(T (z ) = 0, det(T 1 (z, h)) = det(T 2 (z, h)) = 0}. Define varieties W = {(x, y, z, h) | y = y = 0, det(T (z)) = det(T (z )) = 0, x 1 = x 1 , x 2 = x 2 } and W = {(x, y, z, h) | y = y = 0, det(T 1 (z, h)) = det(T 2 (z, h)) = 0}.
Notice that W = W ∩ W . Moreover, we can find a smooth silhouette curve C that is not an orthogonal line to (x 1 , x 2 )-plane and that contains two distinct points q, q , with π C (q) = π C (q ) such that the projection of T q C (resp. T q C) onto R 2 is injective. Notice that j 2

(2) P (q, q ) ∈ W \ W . Hence, W ⊆ W . Moreover, since W is the Cartesian product of determinant varieties (which are irreducible by Proposition 2.1.4(a)) with linear spaces, we have that W is also irreducible [BCR1998, Theorem 2.8.3 (iii)]. In other words, W = W ∩ W is a proper sub-variety of the irreducible variety W . Hence, dim(W ) = dim(W ∩W ) < dim(W ), equivalently, codim(W ) > codim(W ) = 2n. Hence, by Corollary 2.2.16 we have that, generically, the pre-image of W under j 2

(2) P is empty. Since, by Proposition 3.2.12, Assumption A 1 (which is necessary to guarantee that L c is well-defined) is also generic, we imply that, generically, there is no distinct pair q, q ∈ C such that π C (q) = π C (q ) and q ∈ L c , equivalently, L c ∩ L n = ∅ which proves the proposition.

We thus proved the following proposition that the silhouette of a generic surface in R n satisfies all assumptions except for Assumptions A 3(c) and A 5 -that we could not prove. However, based on previous results with three variables [IMP2016b], we formulate the following conjecture.

Proposition 3.2.17. For a generic function

P ∈ C ∞ (R n , R n-2 ), Assumptions A 1 , A 2 , A 3(a) , A 3(b) and A 4 hold.
Conjecture 3.2.18. For a generic function P ∈ C ∞ (R n , R n-2 ), Assumptions A 3(c) and A - 5 hold.

4

Modelling system Our goal in this chapter is, under our weak assumptions (Definition 3.1.1), to encode the singularities of π C (C) by a square and regular (see Definition 4.4.1) system, called the Ball system, so that it is solvable with certified numerical methods.

In Section 4.1, we define this system Ball(P ). Then, we generically study the plane projection of points in L n in Section 4.2. In Section 4.3, we then locally parametrize the curve around the points in L c (Defined in Section 3.1.1) to simplify the computation of Ball(P ) and its Jacobian. In Section 4.4, we determine the necessary and sufficient conditions for this system to be regular.

Encoding the singular points of the plane projection

By our Assumption A - 5 , the singularities of the projected curve π C (C) are only nodes and cusps. Intuitively, a node appears when two points of C project to the same point and a cusp appears when projecting a point with a tangent line orthogonal to the projection plane (see Figure 4.1). The idea to encode the nodes is to design a system whose variables are the coordinates of two different points in R n constrained to be on C and so that they have the same plane projection. To encode a cusp, we design a system whose variables are the coordinates of one point in R n constrained to be on C with a tangent orthogonal to the projection plane.

C is smooth π C (C) Cusps Node (x 1 , x 2 , y) (0, 0, r √ t) y ∈ R n-2 x 1 x 2 (0, 0, -r √ t)
(0, 0, r) Furthermore, in order to apply certified numerical methods we need systems that are square and regular (Definition 4.4.1). To solve this issue and to gather the two systems into a single one, we first parameterize two different points of C with the same projection by (x 1 , x 2 , y + r √ t) and (x 1 , x 2 , y -r √ t), with x 1 , x 2 , t ∈ R, y, r in R n-2 and r = 1, where r denotes the Euclidean norm of r. Then, given any function f from R n to R so that f = 0 is one of the n -1 hypersurfaces that define C, we introduce in Definition 4.1.1 two smooth functions S • f and D • f . When t > 0, they return, roughly speaking, the arithmetic mean and difference of f at the above two points, hence they both vanish if and only if the two points are on the hypersurface f = 0. When t = 0, the two points coincide and S • f and D • f return, roughly speaking, f evaluated at this point and the gradient of f (at that point) scalar the "vertical" vector (0, 0, r); hence, they both vanish if and only if the point is on the hypersurface f = 0 and its tangent hyperplane is normal to the plane of projection. It follows that given a curve defined by P 1 = • • • = P n-1 = 0, the solutions of the so-called Ball system of all S • P i = D • P i = 0 is the set of points on the curve that project to nodes and cusps (Theorem 4.1.3). Note that we consider √ t instead of t in the parameterization (x 1 , x 2 , y ± r √ t) for ensuring the regularity of the Ball system when t = 0 (because this ensures that the linear term of the Taylor expansion of D • f , with respect to t, does not vanish). Definition 4.1.1. Let y, r be two variables in R n-2 and t be a variable in R. For a smooth function f : B ⊂ R n → R, we define the functions:

S •f (x 1 , x 2 , y, r, t) = 1 2 [f (x 1 , x 2 , y + r √ t) + f (x 1 , x 2 , y -r √ t)], for t > 0 f (x 1 , x 2 , y),
for t = 0 and

D•f (x 1 , x 2 , y, r, t) =    1 2 √ t [f (x 1 , x 2 , y + r √ t) -f (x 1 , x 2 , y -r √ t)], for t > 0 ∇f (x 1 , x 2 , y) • (0, 0, r) = n i=3 ∂f ∂x i r i , for t = 0.
Lemma 4.1.2. If f is a smooth function defined on B ⊆ R n , then both S • f and D • f are smooth functions on the subset

B Ball = {(x1, x2, y, r, t) ∈ R × R × R n-2 × R n-2 × R | t 0, (x1, x2, y ± r √ t) ∈ B, r 2 = 1} of R 2n-1
, where r denotes the Euclidean norm of r.

Proof. On the subset B Ball with t > 0, both S • f (x 1 , x 2 , y, r, t) and D • f (x 1 , x 2 , y, r, t) are the compositions of smooth functions, hence they are smooth functions.

For a point X = (x 1 , x 2 , y, r, t) in B Ball with t = 0, we will prove that S

• f (resp. D • f ) is a C s function for an arbitrarily s which implies that S • f (resp. D • f ) is smooth.
First define the function

S 0 • f (x 1 , x 2 , y, r, t) = 1 2 [f (x 1 , x 2 , y + rt) + f (x 1 , x 2 , y -rt)], for t > 0 f (x 1 , x 2 , y),
for t = 0.

Since S 0 • f (x 1 , x 2 , y, r, t) is an even smooth function with respect to t, the partial derivatives of S 0 • f with respect to t of odd orders, evaluated at X, are zero. For an integer s > 0, by the parametrized Taylor formula without remainder [Dem2000, Proposition 4.2.2], there exist smooth functions a i (x 1 , x 2 , y, r),

with integers 0 i < s such that S 0 • f (x 1 , x 2 , y, r, t) = s-1 i=0 a i (x 1 , x 2 , y, r)t 2i + t 2s • φ(x 1 , x 2 , y, t), where φ(x 1 , x 2 , y, t) is a smooth function. Notice that S • f (x 1 , x 2 , y, r, t) = s-1 i=0 a i (x 1 , x 2 , y, r)t i +t s •φ(x 1 , x 2 , y, √ t)
, so that a partial derivative exists up to order s at t = 0. Thus, S • f (x 1 , x 2 , y, r, t) is a C s-1 function. This holds for any arbitrarily large s, hence S

• f (x 1 , x 2 , y, r, t) is a C ∞ function.
Now, we prove that D • f is continuous at X = (x 1 , x 2 , y, r, 0). Let X i be a sequence that converges to X. To prove that D • f (X i ) converges to D • f (X), it is enough to show that for a sequence t i that converges to 0, then we have that D • f (x 1 , x 2 , y, r, t n ) converges to D • f (X). We can assume that t i = 0 for all i, so that

lim t i →0 D • f (x 1 , x 2 , y, r, t i ) = lim t i →0 1 2 √ t i [f (x 1 , x 2 , y + r √ t i ) -f (x 1 , x 2 , y -r √ t i )] = lim t i →0 1 2 √ t i [f (x 1 , x 2 , y + r √ t i ) -(f (x 1 , x 2 , y)- f (x 1 , x 2 , y)) -f (x 1 , x 2 , y -r √ t i )] = lim t i →0 1 2 √ t i [f (x 1 , x 2 , y + r √ t i ) -f (x 1 , x 2 , y)]
+ lim

t i →0 1 2 √ t i [f (x 1 , x 2 , y) -f (x 1 , x 2 , y -r √ t i )] = 1 2 ∇f • (0, 0, r) - 1 2 ∇f • (0, 0, -r) = ∇f • (0, 0, r).
We now prove that D • f is smooth at X. Similarly to the proof of the case of

S • f , since the function 1 2 [f (x 1 , x 2 , y + rt) -f (x 1 , x 2 , y -rt)
] is odd with respect to t, there exist smooth functions b i (x 1 , x 2 , y, r), for 1 i < s and ψ(x 1 , x 2 , y, r, t)

such that 1 2 [f (x 1 , x 2 , y + rt) -f (x 1 , x 2 , y -rt)] = s-1 i=0 b i (x 1 , x 2 , y, r)t 2i+1 + t 2s+1 • ψ(x 1 , x 2 , y, t). Notice that D • f (x 1 , x 2 , y, r, t) = s-1 i=0 b i (x 1 , x 2 , y, r)t i + t s • ψ(x 1 , x 2 , y, √ t)
, so that a partial derivative exists up to order s at t = 0. Thus, D • f (x 1 , x 2 , y, r, t) is a C s-1 function. This holds for any arbitrarily large s, hence

D • f (x 1 , x 2 , y, r, t) is a C ∞ function.
By Lemma 4.1.2, we can build our Ball system based on the smooth operators S and D. We state the main theorem of this section.

Theorem 4.1.3.

Consider P = (P 1 , . . . , P n-1 ) ∈ C ∞ (R n , R n-1 ) satisfy Assump- tions A 1 , A 2 , A 3 and A 4 . Then, X = (x 1 , x 2 , y, r, t) ∈ R × R × R n-2 × R n-2 × R is a solution of the Ball system Ball(P ) =        S • P 1 (X) = • • • = S • P n-1 (X) = 0 D • P 1 (X) = • • • = D • P n-1 (X) = 0 r 2 -1 = 0 (4.1)
if and only if (x 1 , x 2 ) is a singular point of π C (C) (see Definition 4.1.1 for the notation S • P i and D • P i ).

We postpone the proof of Theorem 4.1.3 to the end of Section 4.3. Remark 4.1.4. The purpose behind the last equation r2 -1 = 0 in Ball(P ) is to impose the constraint that r is not trivial. Hence, the equation r 2 -1 = 0 can be replaced with a different constraint. However, this condition is first a polynomial constraint. Second, it provides clear bounds of r when solving Ball(P ). In addition, as a special case, for n = 3, r then consists of one variable. Thus, the equalities 0 = r 2 -1 = r 2 3 -1 imply r 3 = ±1 which reduces the complexity of solving Ball(P ).

We now study a mapping from the solutions of the Ball system to pairs of points on the curve C. Definition 4.1.5. Let P ∈ C ∞ (R n , R n-1 ). Define L n to be the set of pairs (q 1 , q 2 ) with q 1 , q 2 ∈ C, q 1 = q 2 and π C (q 1 ) = π C (q 2 ), also define L c to be the set of pairs (q 1 , q 1 ) with q 1 ∈ L c , and let L = L n ∪ L c . Lemma 4.1.6.

Consider P = (P 1 , . . . , P n-1 ) ∈ C ∞ (R n , R n-1 ) and let X = (x 1 , x 2 , y, r, t) ∈ R × R × R n-2 × R n-2 × R, with r = 1. Assume that P satisfies Assumption A 1 .
Then X is a solution of Ball(P ) if and only if for the points

q 1 = (x 1 , x 2 , y + r √ t) and q 2 = (x 1 , x 2 , y -r √ t), the pair (q 1 , q 2 ) is in L n , or in L c with (0, 0, r) ∈ R × R × R n-2 in T q 1 C.
Proof. Note that, by Assumption A 1 , the tangent space to the curve at any of its points is well defined and is a line. First, assume that X is a solution of Ball(P ). We consider two cases:

(a) If t > 0, then since r = 0 ∈ R n-2 we have that q 1 = q 2 . Moreover, since S •P i (X) = D •P i (X) = 0 for all i ∈ {1, . . . , n-1}, we deduce that P i (q 1 ) = P i (q 2 ) = 0, thus q 1 , q 2 ∈ C. Moreover, since π C (q 1 ) = π C (q 2 ) = (x 1 , x 2 ) we have q 1 , q 2 ∈ L n . Thus, (q 1 , q 2 ) ∈ L n .

(b) If t = 0, then q 1 = q 2 . First, P i (q 1 ) = S • P i (X) = 0, for all indices i ∈ {1, . . . , n -1}, hence q 1 ∈ C. Moreover, we have that 0 = D • P i (X) = ∇P i (q 1 )•(0, 0, r), for all i ∈ {1, . . . , n-1}, equivalently, J P (q 1 )•(0, 0, r) T = 0 ∈ R n-1 , i.e., we have (0, 0, r) ∈ T q 1 C. Thus, q 1 ∈ L c and hence, (q 1 , q 1 ) ∈ L c . Now, let us prove the other direction:

(a) If (q 1 , q 2 ) ∈ L n , then q 1 = q 2 and t = 0. Also, since q 1 , q 2 ∈ C, we can write that S • P i (X) = 1 2 (P i (q 1 ) + P i (q 2 )) = 0, and D

• P i (X) = 1 (b) If (q 1 , q 2 ) ∈ L c and (0, 0, r) ∈ R × R × R n-2 is in T q 1 C, one has q 1 = q 2 ∈ L c ⊆ C
, and t = 0. Moreover, for all i ∈ {1, . . . , n -1} we have S • P i (X) = P i (q 1 ) = 0 and since (0, 0, r) ∈ T q 1 C, we can equivalently write D • P i (X) = ∇P i (q 1 ) • (0, 0, r) = 0. Thus, X is a solution of Ball(P ).

Definition 4.1.7. Let Sol Ball(P ) be the solution set of Ball(P ). Define the function Ω P from Sol Ball(P ) to L that sends X = (x 1 , x 2 , y, r, t) ∈ R×R×R n-2 ×R n-2 ×R to the ordered pair q 1 = (x 1 , x 2 , y + r √ t) and q 2 = (x 1 , x 2 , y -r √ t). Notice that the function Ω P is well-defined by Lemma 4.1.6.

Lemma 4.1.8. If P ∈ C ∞ (R n , R n-1 ) satisfies Assumption A 1 , then Ω P is surjec- tive.
Proof. For any pair (q 1 , q 2 ) ∈ L n we have that the point X = ( 1 2 (q 1 + q 2 ), Π C (q 1 -q 2 ) q 1 -q 2 , 1 4 q 1 -q 2 2 ) ∈ R n × R n-2 × R is a solution of Ball(P ), where Π C (q 1 -q 2 ) is the vector in R n-2 obtained by omitting the first two coordinates (which are zeros) from q 1 -q 2 . Note that Ω P (X) = (q 1 , q 2 ). If the pair (q 1 , q 1 ) is in L c , we define r in the following way, we take a unit vector v ∈ T q 1 C (the first two coordinates of v are zeros since q 1 ∈ L c ). We set r to be Π C (v). Again X = (q 1 , r, 0) ∈ R n × R n-2 × R is a solution of Ball(P ), with Ω P (X) = (q 1 , q 1 ). Thus, Ω P is surjective.

Remark 4.1.9. Notice that if X = (x 1 , x 2 , y, r, t) is in Sol Ball(P ) , then Ω P (X) ∈ L n (resp. Ω P (X) ∈ L c ) if and only if t = 0 (resp. t = 0).

Remark 4.1.10. Preserving the notation in Lemma 4.1.6, notice that if X = (x 1 , x 2 , y, r, t) is a solution of Ball(P ), then X = (x 1 , x 2 , y, -r, t) is another solution. Moreover, both solutions characterize the same unordered pair Ω P (X) = Ω P (X ) = (q 1 , q 2 ). We call X and X twin solutions. An alternative choice would have been to take r in a projective space instead of the sphere to identify these twin solutions.

Example 4.1.11.

Let n = 3 and B = {(x 1 , x 2 , x 3 ) ∈ R 3 | x 1 , x 2 , x 3 ∈ [-2, 2]}. Define P 1 (x 1 , x 2 , x 3 ) = x 1 -(x 3 -1) 3 , P 2 (x 1 , x 2 , x 3 ) = x 2 -(x 3 -1) 2 and P = (P 1 , P 2 ).
The Jacobian matrix of P has full rank over C, thus Assumption A 1 is satisfied. The set L n is empty since π C is injective over C, hence Assumption A 4 is satisfied. The only point of C with a tangent line orthogonal to the (x 1 , x 2 )-plane is q 1 = (0, 0, 1), thus L c = {q 1 } and Assumption A 2 is satisfied. By Lemma 4.3.2, the multiplicity of the system {P = 0, (x 1 , x 2 ) = π C (q 1 )} at its unique solution q 1 is min{ord 1 ((x 3 -1) 3 ), ord 1 ((x 3 -1) 2 )} = min{3, 2} = 2 (ord is defined in Definition 2.2.8). Moreover, for any point q 0 ∈ C different from q 1 , the multiplicity of the corresponding system at its unique solution q 0 is one, thus P satisfies Assumption A 3 . The system Ball(P ):

                  
x 1 -3r 2 ty + 3r 2 t -y 3 + 3y 2 -3y + 1 = 0

x 2 -r 2 t -y 2 + 2y -1 = 0 -r 3 t -3ry 2 + 6ry -3r = 0 -2ry + 2r = 0 r 2 -1 = 0 (4.2)
has two twin solutions X = (0, 0, 1, 1, 0) and X = (0, 0, 1, -1, 0) in B Ball(P ) ⊂ R 2•3-1 = R 5 such that Ω P (X) = Ω P (X ) = (q 1 , q 1 ) ∈ L c .

Example 4.1.12. Let B be defined as in Example 4.1.11. Define the functions P 1 (x 1 , x 2 , x 3 ) = x 1 -(x 2 3 -1), P 2 (x 1 , x 2 , x 3 ) = x 2 -(x 3 3 -x 3 ) and P = (P 1 , P 2 ). The Jacobian matrix of P has full rank over C, thus Assumption A 1 is satisfied. Moreover, the set L c is empty and L n = {q 1 , q 2 }, with q 1 = (0, 0, 1), q 2 = (0, 0, -1), i.e., Assumptions A 2 and A 4 are satisfied. The multiplicity of the system {P = 0 ∈ R n-1 , x 1 = x 2 = 0} at both q 1 , q 2 is equal to one, thus Assumption A 3 is also satisfied. The system Ball(P ):

                   x 1 -r 2 t -y 2 + 1 = 0 x 2 -r 2 ty -y 3 + y = 0 -2ry = 0 -r 3 t -3ry 2 + r = 0 r 2 -1 = 0 (4.3)
has two twin solutions X = (0, 0, 0, 1, 1) and X = (0, 0, 0, -1, 1) in R 5 such that Ω P (X) = Ω P (X ) = (q 1 , q 2 ) ∈ L n .

Singularities induced by L n

Under Assumptions A 1 , A 2 , A 3 and A 4 , we studied in Lemma 3.2.7 the singularity types of the projection of the points in L n . More precisely, for q ∈ L n , the point π C (q) is the intersection of two branches of π C (C), that is, the curve π C (C) at the point π C (q) is either an A 2k+1 singularity with k 0 or equivalent to a curve defined by x 2 1 -g(x 2 2 ) at the origin, with ord(g) = ∞. In both cases, the point π C (q) is a singularity in π C (C). Under our weak or strong assumptions, the singularity type of π C (q) is A 1 (see Definition 2.2.4). We prove in Section 4.4 that this implies the regularity of our Ball system (defined in (4.1)).

Recall that in Lemma 3.2.6 we studied, under Assumptions A 1 , A 2 , A 3 and A 4 , the plane projection of points in C \ (L c ∪ L n ). Hence, in order to complete the classification of the points in π C (C) under the latter assumptions, we study in the next section the plane projection of the points in L c .

Singularities induced by L c

We now study the types of singularities of the plane curve π C (C) obtained by projecting points in L c , that is when the tangent to C is orthogonal to the projection plane. We start by locally parameterizing C around a point in L c . This parameterization will ease the computation of Ball(P ) and its Jacobian in Section 4.4. In the rest of this section and Section 4.4, the analysis is simplified by translating relevant points or assuming the curve C is parameterizable by a specific variable. On the other hand, in the algorithms presented in Chapter 5, the input is not modifiable at all but every computation uses interval arithmetic. This implies that the exact coordinates of a point may not be known, instead we only compute with a box containing it and isolating it from other relevant points. The idea of our semialgorithms is to check that some function does not vanish on such a box. This then implies that such a function does not vanish at the point this box contains. The theoretical analysis of this section can then be applied to the point to deduce the appropriate property without the knowledge of the exact location of that point.

Lemma 4.3.1. Let P ∈ C ∞ (R n , R n-1 ). Let q ∈ L c such that Assumption A 1 is satisfied in a neighbourhood of q in B. Without loss of generality one can assume q = 0 ∈ R n . Then there exist an invertible matrix M of size (n -1) × (n -1) of smooth functions in a neighbourhood of q and smooth functions f 1 , f 2 , f 3 , . . . , f n-1 defined in a neighbourhood of 0 ∈ R, such that:

          x 1 -f 1 (x n ) x 2 -f 2 (x n ) x 3 -f 3 (x n ) . . . . . . x n-1 -f n-1 (x n )           = M •           P 1 P 2 P 3 . . . . . . P n-1           , (4.1) with min{ord(f 1 (x n )), ord(f 2 (x n ))} > 1 (ord is defined in Definition 2.2.8).
Proof. Since rank(J P (q)) = n -1 (Assumption A 1 ), there exists k ∈ {1, . . . , n} such that det(M k (q)) = 0, where M k is the sub-matrix of J P obtained by removing the k-th column. Notice that k ∈ {1, 2}, since q ∈ L c implies that det(M 1 (q)) = det(M 2 (q)) = 0. Without loss of generality, we assume that k = n. Using the implicit function theorem [Corollary 2.7.3][Dem2000], there exist smooth functions f 1 , . . . , f n-1 (see Figure 4.4) of one variable such that we have that

P j (f 1 (x n ), . . . , f n-1 (x n ), x n ) = 0, j ∈ {1, . . . n -1}. (4.2)
Define the function ϕ that maps x i to z i = x i -f i (x n ), for all i ∈ {1, . . . , n -1} and x n to z n = x n . We can see that ϕ is a diffeomorphism and z = (z 1 , . . . , z n ) is a local coordinate system around q. Hence, we can define the function G j (z) = P j • ϕ -1 (z) = P j (x) for all integers 1 j n -1. Using Hadamard's Lemma [Dem2000, Proposition 4.2.3] for the first n -1 variables of z, we can write G j (z) -G j (0, . . . , 0,

z n ) = n-1 i=1 z i • h ji (z) for some smooth func- tions h ji . Note that ϕ -1 (z) = (z 1 + f 1 (z n ), . . . , z n-1 + f n-1 (z n ), z n ). Hence, G j (0, . . . , 0, z n ) = P j • ϕ -1 (0, . . . , 0, z n ) = P j (f 1 (z n ), . . . , f n-1 (z n ), z n ) = P j (f 1 (x n ), . . . , f n-1 (x n ), x n ).
The latter function is equal to zero by (4.2). Thus,

P j (x) = G j (z) = n-1 i=1 z i • h ji (z) = n-1 i=1 (x i -f i (x n )) • H ji (x), with H ji (x) = h ji • ϕ(x). Defining M 0 = H ji 1 j,i n-1 we get:      P 1 . . . . . . P n-1      = M 0 •      x 1 -f 1 (x n ) . . . . . . x n-1 -f n-1 (x n )      .
Notice that M 0 evaluated at q is the invertible matrix M n (q). Hence, by continuity of the determinant function, there is a neighbourhood of q in which M 0 is invertible. Thus, writing M as the inverse of M 0 we get:

Q 0 =      x 1 -f 1 (x n ) . . . . . . x n-1 -f n-1 (x n )      = M •      P 1 . . . . . . P n-1      . (4.3)
To prove that min{ord(f 1 (x n )), ord(f 2 (x n ))} > 1, we take the Jacobian matrices of both sides of (4.3) and we evaluate them at q = 0. We get the equation J Q 0 (q) = M (q) • J P (q). By invertibility of M (q) we deduce that the k-th sub-matrices, obtained by removing the k-th column, of J Q 0 (q) and J P (q) have the same rank.

Computing J Q 0 (q) and considering the fact that det(M 1 (q)) = det(M 2 (q)) = 0 implies that f 1 (0) = f 2 (0) = 0, we thus have that min{ord(f

1 (x n )), ord(f 2 (x n ))} is at least two.
Lemma 4.3.2. Preserving the notation and the assumptions in Lemma 4.3.1, the multiplicity m of the system

S = {Q 0 (x) = 0 ∈ R n-1 , x 1 = x 2 = 0} at q is equal to d = min{ord(f 1 (x n )), ord(f 2 (x n ))}.
Proof. First, we start with the case m < ∞. By Proposition 2.2.11, we can assume, without loss of generality, that f 1 , . . . , f n-1 are polynomials. Following the notation in Definition 2.1.1, let R[x] (resp. R[x n ]) be the ring of polynomials with n variables (resp. one variable) and R[x] q (resp. R[x n ] 0 ) be its localization at q (resp. 0 ∈ R). Also, define I S to be the ideal generated by the polynomials of S in R[x] q (as I G is defined in Definition 2.1.1), i.e.,

I S = x 1 -f 1 (x n ), x 2 -f 2 (x n ), . . . x n-1 -f n-1 (x n ), x 1 , x 2 = x 1 -f 1 (x n ), x 2 - f 2 (x n ), . . . x n-1 -f n-1 (x n ), f 1 (x n ), f 2 (x n ) . If f 1 (x n ) = f 2 (x n ) = 0, then
the ideal I S is of dimension one, hence, S has an infinite number of solutions which contradicts the assumption that m < ∞. Thus, d < ∞ which means that there exist

h 1 , h 2 ∈ R[x n ] 0 such that h 1 (x n )f 1 (x n ) + h 2 (x n )f 2 (x n ) = x d n . Thus, I S = x 1 -f 1 (x n ), x 2 -f 2 (x n ), . . . , x n-1 -f n-1 (x n ), x d n . Note that the set {x 1 -f 1 (x n ), x 2 -f 2 (x n ), . . . x n-1 -f n-1 (x n ), x d n } is a Gröbner basis of I S with respect to Local Lexicographical ordering x 1 > • • • > x n . Hence, by Theorem 2.1.2, we have dim( R[x]q I S ) = dim( R[x]q LT (I S ) ) = dim( R[x]q x 1 ,x 2 ,...x n-1 ,x d n ),
where LT (I S ) is the ideal generated by the leading terms of I S . Consequently,

m = dim( R[x]q I S ) = dim( R[xn] 0 x d n ) = d.
Second, assume that m = ∞. We prove that d = ∞, that is,

∂ k f 1 ∂x k n (0) = ∂ k f 2 ∂x k n
(0) = 0 for any positive integer k. Preserving the notation in Definition 2.2.10, consider

x 3 x 1 x 2 q C (a) The curve C in R 3 where q ∈ L c . x 3 x 1 x1 = f 1 (x3 ) (b)
x 1 is locally a smooth function of x 3 at q.

x 3

x 2

x 2 = f 2 (x 3 ) (c) x 2 is locally a smooth function of x 3 at q. , for c, is non-zero. We consequentially show that

∂ k f 1 ∂x k n (0) = ∂ k f 2 ∂x k n (0) = 0.
We prove the previous statements by induction on k.

For k = 1, since q ∈ L c , we already showed in the proof of Lemma 3.2.3 that a non-trivial element c =

n i=1 v i ∂ ∂x i is in D 1 q [S] \ D 0 q [S] if and only if v = (v 1 , . . . , v n ) is in T q C.
On the other hand, T q C is generated by the vector (f 1 (0), . . . f n-1 (0), 1),

thus c x 1 n = v n = 0. The function f 1 (x n ) is in the set of functions generated by S thus 0 = c • (f 1 (x n )) = n i=1 v i ∂ ∂x i q • (f 1 (x n )) = c x 1 n ∂f 1
∂xn (0), and hence ∂f 1 ∂xn (0) = 0. Thus, the induction hypothesis holds for k = 1.

For k > 1, define c = φ n (c) and consider two cases:

(a) c ∈ D k-1 q [S] \ D k-2 q [S]: By the induction hypothesis, the coefficient c x k-1 n corresponding to ∂ k-1 ∂x k-1 n for c is non-zero and ∂ k f 1 ∂x k n (0) = ∂ k f 2 ∂x k n (0) = 0, for all k < k. Notice that by the definition of φ n , we have c x k n = c x k-1 n = 0. Hence, 0 = c • f 1 (x n ) = k i=1 c x i n ∂ i f 1 ∂x i n (0) = c x k n ∂ k f 1 ∂x k n (0). Hence, ∂ k f 1 ∂x k n (0) = 0. Similarly, we prove that ∂ k f 2 ∂x k n (0) = 0. Thus in Case (a), the lemma is proved. (b) c ∈ D k-2 q [S]: Since c ∈ D k q [S] \ D k-1 q [S], there exists j ∈ {1, . . . , n -1} such that the element c = φ j (c) is in D k-1 q [S] \ D k-2 q [S]. By the induction hypothesis, the coefficient c x k-1 n corresponding to ∂ k-1 ∂x k-1 n for c , is non-zero.
On the other hand,

c x j x k-1 n = c x k-1 n = 0. Hence, since φ n (c x j x k-1 n ∂ k ∂x j ∂x k-1 n ) ∈ D k-1 q [S] \ D k-2 q [S
], then so is φ n (c) = c which contradicts the assumption. Thus, Case (b) is impossible.

With the additional Assumptions A 2 , A 3 and A 4 , one can give a more precise form of f 1 and f 2 in Equation (4.1).

Lemma 4.3.3. Let P ∈ C ∞ (R n , R n-1 ). Let q ∈ L c such that Assumptions A 1 , A 2 , A 3 and A 4 in B, then there exist an invertible matrix M of size (n -1) × (n -1) of smooth functions in a neighbourhood of q, a smooth diffeomorphism ϕ defined in an open subset of R n , with z = (z 1 , . . . , z n ) = ϕ -1 (x) and smooth functions f 3 , . . . , f n-1 , g defined in a neighbourhood of 0 ∈ R, such that

Q =           z 1 -z n • g(z 2 n ) z 2 -z 2 n z 3 -f 3 (z n ) . . . . . . z n-1 -f n-1 (z n )           = M •           P 1 P 2 P 3 . . . . . . P n-1           • ϕ, (4.4)
on a neighbourhood of q. Moreover, either ord(g(z n )) = ∞ or there exists an integer k > 0 with g(z n ) = z k n .

Proof.

Step 1: Equation (4.3) implies that Q 0 and P define the same curve C in a neighbourhood of q and that the function Q 0 satisfies the same assumptions as P around q. By Lemma 4.3.2, d = min{ord(f 1 (x n )), ord(f 2 (x n ))} is the multiplicity of the system {Q 0 (x) = 0 ∈ R n-1 , x 1 = 0, x 2 = 0} at q. By Assumption A 3 , we have that d = 2.

Without loss of generality, assume that ord(f

2 (x n )) = 2 and ∂ 2 f 2 ∂x 2 n (0) = 2. Hence, there is a smooth function v such that f 2 (x n ) = x 2 n (1 + x n • v(x n )). Now, consider the diffeomorphism φ n that sends x n to z n = x n 1 + x n • v(x n ). We have that x 2 -f 2 (x n ) = x 2 -z 2 n . Define f1 (z n ) = f 1 (φ -1 n (z n )) and f2 (z n ) = f 2 (φ -1 n (z n )) = z 2 n . Since ord( f1 (z n )) = ord(f 1 (x n )) d = 2, there exists a smooth function h such that f1 (z n ) = z 2 n h(z n ). Write f1 (z n ) = z 2 n [ h(zn)+h(-zn) 2 + h(zn)-h(-zn)

2

]. Since h(zn)+h(-zn) 2 (resp. h(zn)+h(-zn)

2

) is even (resp. odd), then by Theorem 2.2.2 there exists a smooth function ξ 1 (resp. ξ 2 ) such that h(zn)+h(-zn)

2 = ξ 1 (z 2 n ) (resp. h(zn)-h(-zn) 2 = z n ξ 2 (z 2 n )). Thus, f1 (z n ) = z 2 n (ξ 1 (z 2 n )+z n ξ 2 (z 2 n )).
Notice that ξ 2 (x 2 n ) cannot be the zero function, otherwise f1 ( ) = f1 (-) and f2 ( ) = f2 (-) for all small enough > 0, which contradicts Assumption A 4 .

Step 2: We have two cases:

Case 1: ord(ξ 2 (z n )) = ∞, then define the diffeomorphism φ which sends x 1 to z 1 = x 1 -x 2 ξ 1 (x 2 ), x i to z i = x i for all integers i ∈ {2, . . . , n -1} and x n to z n = x n 1 + x n • v(x n ). Taking g(z n ) = z n ξ 2 (z n ) and ϕ = φ -1 we prove the claim for the first case. Case 2: ord(ξ 2 (z n )) = k < ∞, that is, ξ 2 (z n ) = z k n u(z n ),
for some smooth function u, with u(0) = 0 and an integer k 0. Hence, we can write

x 1 -f1 (z n ) = x 1 -z 2 n ξ 1 (z 2 n ) -z 2k+3 n u(z 2 n ) = x 1 -x 2 ξ 1 (x 2 ) -z 2k+3 n u(x 2 ).
So, defining the diffeomorphism φ which sends x i to z i = x i for all integers i ∈ {2, . . . , n -1},

x n to z n = x n 1 + x n • v(x n ) and x 1 to z 1 = (x 1 - x 2 ξ 1 (x 2 ))u -1 (x 2 ) (which means that x 1 -f 1 (x n ) = u(x 2 )[z 1 -z 2k+3 n ]), we get that:      x 1 -f 1 (x n ) . . . . . . x n-1 -f n-1 (x n )      = u(x 2 ) 0 1×(n-2) 0 (n-2)×1 I n-2 •           z 1 -z 2k+3 n z 2 -z 2 n z 3 -f 3 (z n ) . . . . . . z n-1 -f n-1 (z n )           • φ, (4.5)
for a small enough neighbourhood of q, where I n-2 is the identity matrix of size n -2. Comparing with (4.1), we get:

M •           P 1 P 2 P 3 . . . . . . P n-1           = u(x 2 ) 0 1×(n-2) 0 (n-2)×1 I n-2 •           z 1 -z 2k+3 n z 2 -z 2 n z 3 -f 3 (z n ) . . . . . . z n-1 -f n-1 (z n )           • φ. (4.6) Hence, taking M = u(x 2 ) 0 1×(n-2) 0 (n-2)×1 I n-2 -1
• M and ϕ = φ -1 we recover (4.4).

Following the conclusion of Lemma 4.3.3, the reader may wonder whether the projection of q in π C is always singular. This is clear when g(x n ) = x k n for 0 < k < ∞ since this implies z 2 1 -z k+1 2 = 0 and thus π C (q) is a singularity of the type A 2k . We next prove that the projection is also singular if ord(g(z n )) = ∞. Lemma 4.3.4. Preserving the notation and the assumptions in Lemma 4.3.3, consider the function g defined in (4.4), if ord(g(z

n )) = ∞, then π C (q) is singular in π C (C).
Proof. Since ord(g(z n )) = ∞, then Case 1 in the proof of Lemma 4.3.3 holds. Moreover, we saw in the same proof that ξ 2 (z 2 n ) (restricted to an open neighbourhood of 0 ∈ R) cannot be the zero function. This implies that neither is the function g(z 2 n ) = z 2 n ξ(z 2 n ), i.e., g(z 2 n ), restricted to an open neighbourhood of 0 ∈ R, is not the zero function. Assume for the sake of contradiction that π C (q) is smooth in π C (C), then using the implicit function theorem, there exists a C ∞ -function defined in a neighbourhood of 0 in R, with f (0) = 0 such that for a small neighbourhood of π C (q) in R 2 , one of the following cases is satisfied:

x 1 = x n g(x 2 n ) x 2 = x 2 n ord(g(x n )) = ∞ A 2 k s i n g u l a r i t y , w i t h k > 1 A 2
s i n g u l a r i t y n . Taking the second derivative of both sides with respect to z n and then evaluating at 0 (recall that ord(g(z n )) = ∞), we get the contradiction 0 = 2.

(b) f (z 2 ) = z 1 ⇐⇒ (z 1 , z 2 ) ∈ π C (C). Then f (z 2 n ) = z n g(z 2 n ). The function z n g(z 2
n ) is an odd function but not the zero function, and on the other hand f (z 2 ) is an even function, which leads to a contradiction. Thus, in both cases we have a contradiction, that is, f does not exist and π C (q) cannot be smooth in π C (C).

Returning to (4.4), notice that ϕ is defined in such a way that it preserves the singularity class of π C (C) at the point π C (q). In other words, if C is the plane projection of the curve defined by the Q then (π C (C), 0) and (C, 0) are equivalent. As a corollary of Lemmas 3.2.6, 3.2.7, 4.3.3 and 4.3.4, the points of C in L c ∪ L n are projected to the singular points of π C (C).

Corollary 4.3.5. If P satisfies Assumptions A 1 , A 2 , A 3 and A 4 , then a point q ∈ C projects to a singular point in π C (C) if and only if q ∈ L c ∪ L n .

Proof. If q ∈ L c ∪ L n , then by Lemmas 3.2.7, 4.3.3, and 4.3.4, π C (q) is singular in π C (C). If q ∈ L c ∪ L n , then by Lemma 3.2.6, π C (q) is smooth in π C (C).

Finally, we prove that the solutions of the Ball system project to the singular points of π C (C).

Proof of Theorem 4.1.3: By Corollary 4.3.5, if (x 1 , x 2 ) is singular in π C (C), then there exists a point q 1 ∈ L c ∪L n , with π C (q 1 ) = (x 1 , x 2 ). If q 1 ∈ L c , let q 2 = q 1 and otherwise let q 2 be the unique (by Assumption A 3 ) point in L n , distinct from q 1 , that projects onto (x 1 , x 2 ), i.e. π C (q 1 ) = π C (q 2 ) = (x 1 , x 2 ). Hence, (q 1 , q 2 ) is in L.

Since Ω P is surjective (Lemma 4.1.8), there exists X = (x 1 , x 2 , y, r, t) ∈ Sol Ball(P ) with Ω P (X) = (q 1 , q 2 ).

On the other hand, if X is a solution of Ball(P ), then by Lemma 4.1.6 the pair (q 1 , q 2 ) = Ω P (X) is in L. Hence,

q 1 = (x 1 , x 2 , y + r √ t) ∈ R × R × R n-2 is in L c ∪ L n . Hence, by Corollary 4.3.5 the point (x 1 , x 2 ) is singular in π C (C).

Regularity of the Ball system

In this section, our goal is to prove Theorem 4.4.2 determining, under generic assumptions, necessary and sufficient conditions for Ball(P ) to be regular. We first recall the definition of a regular system. Definition 4.4.1. For some integer m n, let F = (f 1 , . . . , f m ) be a vector of smooth real-valued functions that are defined in R n and let a ∈ R n be a solution of the system {F = 0}. We say that the latter system is regular at a ∈ R n if the rank of its Jacobian matrix, evaluated at a, equals to m. We call {F = 0} regular if it is regular at all of its solutions.

Theorem 4.4.2. Let P ∈ C ∞ (R n , R n-1 ) that satisfies Assumptions A 1 , A 2 , A 3 and A 4 , then P satisfies Assumption A 5 if and only if Ball(P ) is regular in B Ball .

In order to prove Theorem 4.4.2, we are going to show that the Jacobian matrices of Ball(P ) and Ball(Q) evaluated at X have the same rank, where Q is defined in Equation (4.4). Recall that Equation (4.4) implies that P and Q define the same curve around q. Notice also that if X = (q, r, 0) ∈ R n ×R n-2 ×R is in Ω -1 P ((q, q)), then X ∈ Ω -1 Q ((q, q)).

Lemma 4.4.3. Let P and Q be as defined in (4.4). Under Assumption A 1 , let (q, r, 0) ∈ R n ×R n-2 ×R be a solution of the system Ball(P ) in B Ball , then Ball(P ) is regular at (q, r, 0) if and only if Ball(Q) is regular at the point (0, r, 0) ∈ R n × R n-2 × R (recall that for simplicity, we assume in Lemma 4.3.3 that q = 0 ∈ R n ).

Proof. Let us write X = (q, r, 0). We are going to prove that the Jacobian matrices of Ball(P ) and Ball(Q) evaluated at X have the same rank. By Remark 4.1.9 we have that Ω P (X) = (q, q) ∈ L c (see Definitions 4.1.7 and 4.1.5), and hence, q ∈ L c . By Lemma 4.1.6 we have that (0, 0, r) ∈ T q C. We prove the claim in three steps:

Step 1: Let M = (f ij ) 1 i,j n-1 be as defined in the Equality (4.4). We define S • M (resp. D • M ) to be the matrix

(S • f ij ) 1 i,j n-1 (resp. (D • f ij ) 1 i,j n-1 ). Using the identity 1 2 (ab + cd) = 1 4 (a + c)(b + d) + 1 4 (a -c)(b -d), one deduces the properties for any f, g ∈ C ∞ (R n , R): S • f g = (S • f )(S • g) + t(D • f )(D • g) (4.1) D • f g = (D • f )(S • g) + (S • f )(D • g) (4.2)
These identities applied to Equation (4.4) yield

     S • Q 1 . . . . . . S • Q n-1      = S • M tD • M •               S • (P 1 • ϕ) . . . . . . S • (P n-1 • ϕ) D • (P 1 • ϕ) . . . . . . D • (P n-1 • ϕ)               and      D • Q 1 . . . . . . D • Q n-1      = D • M S • M •               S • (P 1 • ϕ) . . . . . . S • (P n-1 • ϕ) D • (P 1 • ϕ) . . . . . . D • (P n-1 • ϕ)              
Combining the last two equalities:

              S • Q 1 . . . . . . S • Q n-1 D • Q 1 . . . . . . D • Q n-1               = S • M tD • M D • M S • M •               S • (P 1 • ϕ) . . . . . . S • (P n-1 • ϕ) D • (P 1 • ϕ) . . . . . . D • (P n-1 • ϕ)               (4.3) Notice that S • M tD • M D • M S • M X = M (q) 0 D • M (X) M (q)
(recall that in our case we have S • M (X) = M (q) ) and that the latter matrix has an inverse (recall that, by Lemma 4.3.3, M (q) is an invertible matrix of size n -1), namely,

M (q) -1 0 -M (q) -1 • (D • M )(X) • M (q) -1 M (q) -1 which implies (by continuity of the determinant function) that S • M tD • M D • M S • M is invertible in a neighbourhood of X.
Step 2: Writing y = (y 3 , . . . , y n ) and r = (r 3 , . . . , r n ), consider the diffeomorphism ϕ defined in Lemma 4.3.3 and define the smooth function ψ over an open subset of R 2n-1 containing X which maps the point (x 1 , x 2 , y, r, t) to (ϕ 1 , ϕ 2 , S • ϕ 3 , . . . , S • ϕ n , D • ϕ 3 , . . . , D • ϕ n , t). Notice that we have:

S • (P j • ϕ) = (S • P ) • ψ and D • (P j • ϕ) = (D • P ) • ψ, for 1 j n -1, (4.4) since ϕ i (x 1 , x 2 , y ± r √ t) = ψ i ± ψ n+i-2 ψ 2n-1
for all i ∈ {3, . . . , n}. In fact, using the last two equalities we can also see that ψ -1 exists and is smooth. Thus, ψ is a diffeomorphism.

Step 3: Now, comparing (4.3) with (4.4) we get:

SD • Q :=               S • Q 1 . . . . . . S • Q n-1 D • Q 1 . . . . . . D • Q n-1               = S • M tD • M D • M S • M •               S • P 1 . . . . . . S • P n-1 D • P 1 . . . . . . D • P n-1               • ψ.
Consider the vector SD • P = (S • P 1 , . . . , S • P n-1 , D • P 1 , . . . , D • P n-1 ) T and let J SD•P , J SD•Q and J ψ be the Jacobian matrices of SD •P, SD •Q and ψ respectively. Taking the Jacobian matrix of both sides of the last equality yields:

J SD•Q = S • M tD • M D • M S • M •J SD•P •J ψ +Jacobian( S • M tD • M D • M S • M )•               S • P 1 . . . . . . S • P n-1 D • P 1 . . . . . . D • P n-1               •ψ.
Evaluating the last equality at X = (0, r, 0) and using the fact that ψ(X) = ψ(0, r, 0) = (0, r, 0) = X, we note that the second term of the right-hand side is zero. One thus has:

J SD•Q (X) = S • M tD • M D • M S • M X • J SD•P (X) • J ψ (X). (4.5) Computing J ψ (X), we get J ψ (X) = ∂ϕ 1 ∂z 1 (0) ∂ϕ 1 ∂z 2 (0) 0 1×(2n-3) 0 (2n-2)×1 I 2n-2 , with ∂ϕ 1
∂z 1 (0) = 0 according to the formula in Lemma 4.3.3. Hence by Equation (4.5), it is straightforward to check that:

J Ball(Q) = J SD•Q (X) 2X =    S • M tD • M 0 D • M S • M 0 0 1×(n-1) 0 1×(n-1) 1    X • J SD•P (X) 2X • J ψ (X) =    S • M tD • M 0 D • M S • M 0 0 1×(n-1) 0 1×(n-1) 1    X • J Ball(P ) (X) • J ψ (X). (4.6) Recalling that J ψ (X) and    S • M tD • M 0 D • M S • M 0 0 1×(n-1) 0 1×(n-1) 1    X
are invertible matrices, the proof of the lemma follows. Now, we are ready to prove Theorem 4.4.2 which characterizes the regularity of the solutions of Ball(P ) under generic assumptions. We split the proof in two Lemmas 4.4.6 and 4.4.7. Before that, we introduce a new assumption that helps to simplify the proof. Definition 4.4.4. Let (q 1 , q 2 ) ∈ L. We say that (q 1 , q 2 ) satisfies Assumption A - 5 if q 1 and q 2 are isolated in L n ∪ L c and any of the following conditions is satisfied: (a) If (q 1 , q 2 ) ∈ L n , then the plane projections of the tangent lines of q 1 and q 2 to C are linearly independent.

(b) If (q 1 , q 2 ) ∈ L c , then the plane projection of a small enough neighbourhood of q 1 in C is an ordinary cusp at π C (q 1 ) and the multiplicity of the system {P (x) = 0, (x 1 , x 2 ) = π C (q 1 )} at q 1 is two.

Remark 4.4.5. Assumption A - 5 can be seen as a "local version" of Assumption A - 5 . We are going to prove that if Assumptions A 1 , A 2 , A 3 and A 4 are satisfied, then Assumption A - 5 is equivalent to the condition that Assumption A - 5 is satisfied for all L.

The main reason behind introducing Assumption A - 5 , is that we are going to prove in Lemma 4.4.6 that, under Assumption A 1 , a pair (q 1 , q 2 ) ∈ L satisfies Assumption A - 5 if and only if every X in Ω -1 P ((q 1 , q 2 )) is a regular solution of Ball(P ), whereas Assumption A - 5 is, in general, not sufficient for the regularity of the solutions of Ball(P ). For example, take n = 3 and P = (x 1 -(x 3 -1) 6 , x 2 -(x 3 -1) 9 ) (Figure 4.6). We can see that P satisfies Assumption A 1 , the set L c consists of a unique point q = (0, 0, 1) and the set L n is empty. The plane projection of C is the curve given by the equation x 3 1 -x 2 2 = 0. Hence, Assumption A - 5 is satisfied. However, the multiplicity of the system S = {P (x 1 , x 2 , x 3 ) = 0 ∈ R 2 , x 1 = x 2 = 0} at the point q equals 6 (Lemma 4.3.2). Hence, Assumption A - 5 is not satisfied and one can also check that Ball(P ) is not regular. Lemma 4.4.6. Let P ∈ C ∞ (R n , R n-1 ) that satisfies Assumption A 1 . Let X be a solution of Ball(P ) and (q 1 , q 2 ) = Ω P (X) (Definition 4.1.7), then X is a regular solution of Ball(P ) if and only if (q 1 , q 2 ) satisfies Assumption A - 5 .

Proof.

Let X = (x 1 , x 2 , y, r, t) ∈ R × R × R n-2 × R n-2 ×
R be a solution of Ball(P ). We consider two cases:

Case (a): t = 0, i.e., q 1 = q 2 . It is easy to see that

∂(S•P i ) ∂x j , ∂(D•P i ) ∂x j , ∂(S•P i ) ∂r k , ∂(D•P i ) ∂r k , ∂(S•P i ) ∂t , ∂(D•P i ) ∂t are, respectively, equal to: S • ∂(P i ) ∂x j , D • ∂(P i ) ∂x j , t • D • ∂(P i ) ∂x k , S • ∂(P i ) ∂x k , 1 2 n m=3 D • ( ∂P i ∂xm ) • r m , 1 2t [ n m=3 S • ( ∂P i ∂xm ) • r m -D • P i ].
Hence, by computing the Jacobian matrix of the Ball(P ) we get the matrix:

C C π C(C ) = π C (C ) q x 1 x 3 Figure 4
.6: Comparing the curve C defined in Remark 4.4.5 with C that is defined by x 1 -(x 3 -1) 2 = x 2 -(x 3 -1) 3 = 0 where both of them have the same plane projection which contains an ordinary cusp; however, the multiplicities of S (defined in the same remark) at q differ. In this example, the pair (q, q) satisfies Assumption A - 5 for C , however, it does not for C.

                         S • ∂(P 1 ) ∂x 1 . . . S • ∂P 1 ∂xn t • D • ∂(P 1 ) ∂x 3 . . . t • D • ∂(P 1 ) ∂xn 1 2 n m=3 D • ( ∂P 1 ∂xm ) • rm • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • S • ∂(P n-1 ) ∂x 1 . . . S • ∂(P n-1 ) ∂xn t • D • ∂(P n-1 ) ∂x 3 . . . t • D • ∂(P n-1 ) ∂xn 1 2 n m=3 D • ( ∂P n-1 ∂xm ) • rm D • ∂(P 1 ) ∂x 1 . . . D • ∂(P 1 ) ∂xn S • ∂(P 1 ) ∂x 3 . . . S • ∂(P 1 ) ∂xn 1 2t [ n m=3 S • ( ∂P 1 ∂xm ) • rm -D • P 1 ] • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • D • ∂(P n-1 ) ∂x 1 . . . D • ∂(P n-1 ) ∂xn S • ∂(P n-1 ) ∂x 3 . . . S • ∂(P n-1 ) ∂xn 1 2t [ n m=3 S • ( ∂P n-1 ∂xm ) • rm -D • P n-1 ] 0 . . . 0 2r 3 . . . 2rn 0                          .
We denote by C i (resp. L i ) the i-th column (resp. line) of the latter matrix. Replace

the last column C 2n-1 with n-2 m=1 r m+2 2t C n+m + C 2n-1
, also for all integers 1 k n -1 we replace the line L k with L k + √ t • L k+n-1 and then the line L k+n-1 with L k -2 √ tL k+n-1 . The resulting matrix is:

                 ∂(P 1 )
∂x 1 (q 1 ) . . . ∂P 1 ∂xn (q 1 ) (t) • ∂(P 1 ) ∂x 3 (q 1 ) . . . (t) ∂(P 1 ) ∂xn (q 1 ) 0

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ∂(P n-1 )
∂x 1 (q 1 ) . . . ∂P n-1 ∂xn (q 1 ) (t) • ∂(P n-1 ) ∂x 3 (q 1 ) . . . (t) ∂(P n-1 ) ∂xn (q 1 )) 0

∂(P 1 )
∂x 1 (q 2 ) . . . ∂(P 1 ) ∂xn (q 2 ) -(t) ∂(P 1 ) ∂x 3 (q 2 ) . . .

-(t) ∂(P 1 ) ∂xn (q 2 ) 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ∂(P n-1 ) ∂x 1 (q 2 ) . . . ∂(P n-1 ) ∂xn (q 2 ) -(t) ∂(P n-1 ) ∂x 3 (q 2 ) . . . -(t) ∂(P n-1 ) ∂xn (q 2 )) 0 0 . . . 0 2r 3 . . . 2r n 1 2t                 
.

The determinant of the latter matrix is zero if and only if the determinant of the following matrix is zero:

M 0 = N P (q 1 ) M P (q 1 ) M P (q 1 ) N P (q 2 ) M P (q 2 ) -M P (q 2 ) , where M P (q 1 ), M P (q 2 ) are the sub-matrices that are obtained, respectively, by removing the first two columns from J P (q 1 ), J P (q 2 ) and N P (q 1 ), N P (q 2 ) are the matrices formed by the first two columns of J P (q 1 ), J P (q 2 ), respectively. By linear operations on M 0 , we can see that M 0 has same rank as the matrix M (q 1 , q 2 ) (see Definition 3.2.8). Thus, X is regular for Ball(P ) if and only if M (q 1 , q 2 ) is invertible. By Lemma 3.2.9 we have that M (q 1 , q 2 ) is invertible if and only if none of q 1 , q 2 is in L c (and hence none of the plane projections of T q 1 C, T q 2 C is trivial) and the plane projection of their tangent spaces are different. Equivalently, the pair (q 1 , q 2 ) is in L n and satisfies Assumption A - 5 . Case (b): t = 0, i.e., q 1 = q 2 . Let us write q = q 1 . We prove the claim in three steps:

Step 1: We first simplify P . Without loss of generality and by Lemma 4.3.1we can assume that q = 0 and P 1 , . . . , P n-1 are, respectively, equal to

x 1 -f 1 (x n ), x 2 - f 2 (x n ), . . . , x n-1 -f n-1 (x n ) with the property that min{mult(f 1 ), mult(f 2 )} 2.
For all i ∈ {3, . . . , n -1}, using Taylor's theorem, we can write

f i (x n ) = 3 j=1 a i,j x j n + x 4 n h i (x n ), for some a i,j ∈ R and smooth functions h i (x n ). Since min{mult(f 1 ), mult(f 2 )} 2, we can write f 1 (x n ) = 3 j=2 α j x j n + x 4 n h 1 (x n ) and f 2 (x n ) = 3 j=2 β j x j n + x 4 n h 2 (x n ). Notice that (f 1 (x n ), f 2 (x n ), f 3 (x n ), . . . , f n-1 (x n ), x n )
is a local parameterization system of C around q. Since dim(T q C) = 1 (Assump-tion A 1 ), there exists λ ∈ R * with (a 3,1 , . . . , a n-1,1 , 1) = λr (because the vectors (0, 0, r) ∈ R × R × R n-2 and (0, 0, a 1,3 , . . . , a 1,n-1 , 1) are in T q C \ {0}). In particular, r n = 0.

Step 2: Now, we compute J Ball(P ) at X = (x 1 , x 2 , y, r, 0) by first computing it for X t = (x 1 , x 2 , y, r, t) with t = 0, and then taking the limit when t goes to 0. Since the operator S is linear, we write S(

x i -f i (x n )) = S(x i - 3 j=1 a i,j x j n )-S(x 4 n h i (x n )).
On the other hand, using the identity (4.1) we deduce that S(x 4

n h i (x n )) = S(x 4 n ) • S(h i (x n )) + tD(x 4 n ) • D(h i (x n )), for all ∈ {1, . . . , n -1}. It is straightforward to see that S(x 4 n ) = r 4 n t 2 + 6r 2 n tx 2 n + x 4
n and tD(x 4 n ) = 4r 3 n x n t 2 + 4r n x 3 n t with r = (r 3 , . . . , r n ). Hence, all of the first-order partial derivatives of S(x 4 n h i (x n )), evaluated at X t , converge to zero when t goes to 0. Hence, the partial derivatives of the functions S(x i -f i (x n )) and S(x i -3 j=1 a i,j x j n ) evaluated at X are equal. Using an analogous argument, we deduce that the evaluation of the partial derivatives of the functions D(x i -f i (x n )) and D(x i -3 j=1 a i,j x j n ), at X are also equal. Thus, J Ball(P ) (X t ) and J Ball(P ) (X t ) converge to the same limit J Ball(P ) (X), where P is the function obtained by truncating P beyond degree 3 with respect to the variable x n .

Computing J Ball(P ) (X) = lim t→0 J Ball(P ) (X t ), we get:

                          1 0 . . . 0 0 0 . . . 0 -α 2 r 2 n 0 1 . . . 0 0 0 . . . 0 -β 2 r 2 n 0 0 . . . . . . -a 3,1 0 . . . 0 -a 3,2 r 2 n • • • • • • . . . • • • • • • • • • • • • • • • . . . • • • • • • • • • 0 0 . . . 1 -a n-1,1 0 . . . 0 -a n-1,2 r 2 n 0 0 . . . . . . -2α 2 r n 0 . . . 0 -α 3 r 3 n 0 0 . . . . . . -2β 2 r n 0 . . . 0 -β 3 r 3 n 0 0 . . . . . . -2a 3,2 r n 1 . . . -a 3,1 -a 3,3 r 3 n • • • • • • . . . • • • • • • • • • • • • • • • • • • . . . • • • • • • 0 0 . . . . . . -2a n-1,2 r n 0 . . . 1 -a n-1,1 -a n-1,3 r 3 n 0 0 . . . . . . 0 2r 3 . . . 2r n-1 2r n 0                           .
Hence, observing that the matrix is block diagonal, its determinant is zero if and only if the determinant of the following one is:

          -2α 2 r n 0 . . . 0 0 -α 3 r 3 n -2β 2 r n 0 . . . 0 0 -β 3 r 3 n -2a 3,2 r n 1 0 . . . 0 -a 3,1 -a 3,3 r 3 n . . . . . . -2a n-1,2 r n 0 0 . . . 1 -a n-1,1 -a n-1,3 r 3 n 0 2r 3 . . . 2r n-1 2r n 0          
.

Shifting the columns of the last matrix we get:

          -α 3 r 3 n -2α 2 r n 0 . . . 0 0 -β 3 r 3 n -2β 2 r n 0 . . . 0 0 -a 3,3 r 3 n -2a 3,2 r n 1 0 . . . 0 -a 3,1 . . . . . . -a n-1,3 r 3 n -2a n-1,2 r n 0 0 . . . 1 -a n-1,1 0 0 2r 3 . . . 2r n-1 2r n           .
To compute the determinant of the second block, we expand it about the last row. Hence, the determinant of the last matrix is zero if and only if

r n (α 2 β 3 -α 3 β 2 )(r n + n-1 i=3 a i,1 r i ) = 0. Notice that, by
Step 1, we have that r n = 0 and the third factor

(r n + n-1 i=3 a i,1 r i ) is never zero since it is equal to λ. Thus, J Ball(P ) (X) is invertible iff α 2 β 3 -α 3 β 2 = 0, equivalently, the matrix A = α 2 α 3 β 2 β 3 is invertible.
Step 3: We now show that the invertibility of A is equivalent to the condition that (q, q) satisfies Assumption A - 5 . First assume that A is invertible. It follows that either α 2 = 0 or β 2 = 0 and this yields that the minimum of the multiplicities of f 1 and f 2 is 2. By Lemma 4.3.2, the multiplicity of the system

{P (x 1 , x 2 , y) = 0 ∈ R n-1 , (x 1 , x 2 ) = π C (q)} at q is equal to 2, thus Assumption A - 5 (b) is satisfied. Using the same notation as in the proof of Lemma 4.3.3, one can write f1 (z n ) = z 2 n (ξ 1 (z 2 n ) + z n ξ 2 (z 2 n )). Notice that ξ 2 (x 2
n ) cannot be the zero function, otherwise f1 ( ) = f1 (-) and f2 ( ) = f2 (-) for all small enough > 0, which means that X would be the limit of solutions X of Ball(P ) with Ω P (X ) ∈ L n . X would then be a non-isolated solution and thus a non-regular solution of Ball(P ) which contradicts the assumption. We then have two cases as in Lemma 4.3.3. The first one is when mult(ξ 2 (z n )) = ∞, that would imply that α 2 = α 3 = 0 and contradicts the invertibility of A. We then must satisfy the second case mult(ξ 2 (z n )) = k < ∞ and, after a change of variables, the first equation of the system becomes equivalent to z 1 -z 2k+3 n = 0. The invertibility of A implies that k = 0. The projection of the curve in the plane is thus locally parameterized by (z 3 n , z 2 n ) and is an ordinary cusp, Assumption A - 5 (a) is satisfied. Second, assume that Assumption A - 5 is satisfied. By Lemma 4.3.2 and Assumption A - 5 (b), the minimum of the multiplicities of f 1 and f 2 is 2. Using the proof of Lemma 4.3.3 once again, one can assume that

f 2 (z n ) = z 2 n and f 1 (z n ) = z n g(z 2 n ) or f 1 (z n ) = z 2k+3 n
. By Assumption A - 5 (a), the projection is an ordinary cusp and thus has a parameterization of the form

(z 2 n , z 3 n ), that is f 1 (z n ) = z 3 n .
(q 1 , q 2 ) ∈ L ⊂ B × B.
Proof. Assume that Assumption A - 5 is satisfied and (q 1 , q 2 ) ∈ L. If (q 1 , q 2 ) ∈ L c , then by Lemma 4.3.3 and Assumption A - 5 we must have that the plane projection of a small enough neighborhood of q 1 in C is an ordinary cusp at π C (q 1 ). By Assumption A 3 and Lemma 3.2.3, the multiplicity of the mentioned system at q 1 = q 2 is two. Thus, (q 1 , q 2 ) satisfies Assumption A - 5 . If (q 1 , q 2 ) ∈ L n , then by Lemma 3.2.7 and Assumption A - 5 , we have that π C (q 1 ) is a node in π C (C). Thus, we have that π C (q 1 ) is a transverse intersection of two smooth branches of π C (C). Those branches are the plane projections of two disjoint branches of C each of which contains either q 1 or q 2 . Hence, the plane projections of the tangent spaces of q 1 and q 2 to C are linearly independent. Thus, (q 1 , q 2 ) satisfies Assumption A - 5 . Assume conversely that A - 5 is satisfied for all (q 1 , q 2 ) ∈ L. By Corollary 4.3.5, any singular point of π C (C) is the plane projection of a point q 1 ∈ L c ∪ L n . For some q 2 ∈ C, the pair (q 1 , q 2 ) is in L (which satisfies Assumption A - 5 ). Hence, if (q 1 , q 2 ) is in L n (resp. in L c ) the plane projection of q 1 is a node (resp. an ordinary cusp) by Lemma 3.2.7 (resp. Lemma 4.3.3). Lemmas 4.4.6 and 4.4.7 then imply Theorem 4.4.2.

In Chapter 7, we present examples that show how the Ball system encodes the singularities of π C (C). The next chapter is dedicated to provide a method to examine whether a given curve satisfies the weak or strong assumptions.

5

Semi-algorithms to check the assumptions and isolate singularities In this chapter we present Semi-algorithm 3 that checks the weak assumptions of Section 3.1.2 and, if it terminates, outputs a superset of isolating boxes of the singularities of π C (C). We also present Semi-algorithm 4 that checks the strong assumptions of Section 3.1.2 and, if it terminates, outputs a set of isolating boxes of the singularities of π C (C). The main idea of these semi-algorithms comes from Theorems 4.1.3 and 4.4.2: Theorem 4.1.3 states that, under Assumptions A 1-4 , the singularities of π C (C) are the plane projections of the solutions of Ball(P ). Theorem 4.4.2 states that, under the further Assumption A - 5 , Ball(P ) is regular, so we can use certified numerical methods such as interval Newton methods [MKC2009] to solve Ball(P ). In addition, in order to verify these assumptions, we use subdivision approaches based on interval arithmetic in a semi-algorithm framework.

These semi-algorithms are theoretically relevant but they are not very practical. We thus present in Section 5.3 practical improvements and Semi-algorithm 6.

We present in Section 5.1 the basics of interval arithmetic with the notation and definitions by Lin and Yap [LY2011] and Semi-algorithm 3 and 4 in Section 5.2.

Interval arithmetic

Recall that for some positive integer k, by a closed (resp. open) k-box B, we mean the Cartesian product of k closed (resp. open) intervals. The width of a box B, denoted by w(B), is the maximal length of the intervals of that product. For a subset A ⊂ R k , the set IA is the set of all closed k-boxes that are contained in A.

Functions over boxes. ([MKC2009, §2.3]) closed n-boxes can be seen as an extension of R n . Hence, it would be interesting to speak about extending a function f defined on R n to a function f over to IR n , preserving the "behavior" of a f . More formally: Definition 5.1.1. For the positive integer m, A ⊆ R k , and a function f :

A → R m , the function f : IA → IR m is called an inclusion of f if the set f (B) = {f (x) | x ∈ B} is contained in f (B), for all B ∈ IA. An inclusion f of f is called a box function, if for any descending sequence of closed k-boxes B 1 ⊇ B 2 ⊇ . . . that converges to a point q ∈ R k (that is, i B i = {q}), the sequence f (B 1 ) ⊇ f (B 2 ) ⊇ . . . converges to f (q).
Example 5.1.2. For the case of intervals, the first examples that come to mind are inclusions of arithmetic operations. For instance, the addition +, can be seen as a function over R 2 to R. Hence, the addition of two intervals B 1 = [a, b] and B 2 = [c, d] can be defined to be the interval

B 1 + B 2 = [a + c, b + d].
This inclusion is a box function. The same can be said about multiplication and subtraction. Of course, one must be careful when defining the division since the dominator has to be zero-free. Addition and subtraction can be generalized over higher-dimensional boxes in a natural way. Moreover, outward rounding of an interval ([MKC2009, §3.2]) means to move the left endpoint to the left and the right endpoint to the right at the last digit carried. Outward rounding, when evaluating a box function, guarantees that these enclosures are rigorous when using fixed-point arithmetic (representation of real numbers approximation with a fixed number of digits). We are going to assume that outward rounding is used for any box evaluation performed.

In the rest of this chapter, we assume that we are given a box function f for any function f we consider. The command subdivide is applied to a closed k-box B, and it returns the set of boxes obtained by bisecting B in all dimensions.

An interval matrix M is a matrix whose coefficients are intervals. It can also be seen as the set of all matrices whose (i, j)-th coefficients belong to the (i, j)th interval. The rank of an interval matrix M , denoted by rank( M ), is the minimum of the ranks of all the matrices in this set.

Semi-algorithms to check assumptions and isolate singularities

This section is dedicated to prove the following theorem. Recall that the weak and strong assumptions are defined in Definition 3.1.1.

Theorem 5.2.1. For an open n-box B and a smooth function P from B to R n-1 , Semi-algorithm 3 stops if and only if P satisfies the weak assumptions in B and then it returns a set of isolating boxes of all the singularities of π C (C), plus possibly other spurious disjoint boxes. Semi-algorithm 4 stops if and only if P satisfies the strong assumptions in B and then it returns a set of isolating boxes of all the singularities of π C (C).

To check whether a given function P satisfies the weak assumptions (A 1 , A 2 , A 3 , A 4 and A - 5 ) in B, we use their relation to the solutions of Ball(P ) studied in the previous sections. Recall that for any subset A ⊆ R n , we defined

A Ball = {(x 1 , x 2 , y, r, t) | t 0, (x 1 , x 2 , y + r √ t), (x 1 , x 2 , y -r √ t) ∈ A, r 2 = 1}. Let B be
an open n-box and P be a smooth function from B to R n-1 that satisfies Assumption A 1 in B. Consider the following assumptions: ℵ 1 -All solutions of Ball(P ) in B Ball are regular.

ℵ 2 -For every solution X of Ball(P ) in B Ball , none of the points of the pair Ω P (X) (Definition 4.1.7) is in the boundary of B.

ℵ 3 -No two distinct solutions of Ball(P ) in B Ball , except the twin solutions (Remark 4.1.10), have the same plane projection.

The next lemma shows the relation between these new assumptions and those of Section 3.1.2. The motivation of these alternative assumptions is that they are stated in terms of Ball(P ), which makes them easier to verify in our semi-algorithms.

Lemma 5.2.2. Let B be an open n-box and P be a smooth function from B to R n-1 that satisfies Assumption A 1 in B. Then, Assumptions ℵ 1 , ℵ 2 and ℵ 3 are satisfied if and only if Assumptions A 2 , A 3 , A 4 and A - 5 are satisfied in B.

Proof. If Assumptions A 2 , A 3 , A 4 and A - 5 are satisfied in B, then by Theorem 4.4.2 we have Assumption ℵ 1 is satisfied. Moreover, by Assumptions A 2 and A 4 we have that none of L n , L c intersects ∂B. By Definition 4.1.7, for any solution X of Ball(P ), we have that the points of the pair Ω P (X) are in L n ∪ L c and hence are not in ∂B which implies that Assumption ℵ 2 is satisfied. Assume that Assumption ℵ 3 is not satisfied, that is, there exist two distinct non-twin solutions X, X that have the same plane projection p ∈ R 2 . Let (q 1 , q 2 ) = Ω P (X) and (q 1 , q 2 ) = Ω P (X ). By Lemma 4.1.6, the pairs (q 1 , q 2 ), (q 1 , q 2 ) are distinct and the points q 1 , q 2 , q 1 , q 2 have the same plane projection p. By Assumption A 3 , we cannot have three pairwise distinct points among q 1 , q 2 , q 1 , q 2 . Moreover, if the multiplicity at all of the points q 1 , q 2 , q 1 , q 2 is one, then (q 1 , q 2 ), (q 1 , q 2 ) are in L n and not distinct. Hence, at least a point say q 1 has multiplicity larger than one, i.e., q 1 ∈ L c (Lemma 3.2.3). Hence, the number of solutions counted with multiplicity is at least three which contradicts Assumption A 3 . Hence, Assumption ℵ 3 is satisfied. Now, assume that Assumptions ℵ 1 , ℵ 2 and ℵ 3 are satisfied. Since, by Assumption ℵ 1 , Ball(P ) is a regular square system, its solution set is a zero-dimensional manifold in the compact set B Ball(P ) (regular value theorem). Hence, Ball(P ) has a finite number of solutions in B Ball . Since Ω P (Definition 4.1.7) is surjective (Lemma 4.1.8), the set L (Definition 4.1.5) is also finite. Hence, the set L c ∪ L n is finite (since L c ∪ L n is the image of L under the surjective function (q 1 , q 2 ) → q 1 ). Moreover, by Assumption ℵ 2 , the set L n ∪L c does not intersects the boundary of B. Hence, Assumptions A 2 and A 4 are satisfied in B. To prove that Assumption A 3 is satisfied, let p = (α, β) ∈ π C (C) and |π -1 (p)| 3. For pairwise distinct points q 1 , q 2 , q 3 ∈ π -1 (p), by Lemma 4.1.6, we have that there exist two distinct non-twin solutions X, X of Ball(P ), with Ω P (X) = (q 1 , q 2 ) and Ω P (X ) = (q 1 , q 3 ) such that we have the same plane projection p which contradicts Assumption ℵ 3 . Hence, π -1 C (p) consists of at most two distinct points. We consider two cases:

(a) π -1 C (p) has two distinct elements, say q 1 , q 2 . By Lemma 4.1.6, the pair (q 1 , q 2 ) is in L n , and hence, there exists a solution X = (α, β, y, r, t) ∈ R × R × R n-2 × R n-2 × R of Ball(P ), with t = 0 and Ω P (X) = (q 1 , q 2 ). Since X is a regular solution (Assumption ℵ 1 ), by Lemma 4.4.6 we have that none of q 1 , q 2 is in L c . Hence, by Lemma 3.2.3, the multiplicity of {P

(x 1 , x 2 , y) = 0 ∈ R n-1 , x 1 -α = x 2 -β = 0} at q 1 (resp. q 2 )
is one. Thus, the number of solutions counted with multiplicity is two.

(b) π -1
C (p) has a unique point q. Let m denote the multiplicity of the system

{P (x 1 , x 2 , y) = 0 ∈ R n-1 , x 1 -α = x 2 -β = 0} at q. If m = 1, then
we are done. If m > 1, then by Lemma 3.2.3 we have that q ∈ L c . Hence, there exists a solution of Ball(P ) of the form X = (α, β, y, r, 0) ∈ R × R × R n-2 ×R n-2 ×R such that Ω P (X) = (q, q) (Lemma 4.1.8). Since X is regular (Assumption ℵ 1 ), by Lemma 4.4.6 we have that (q, q) satisfies assumption A 5 .

In particular, the multiplicity m is equal to two.

Thus, for all p ∈ π C (C) the sum of the multiplicities of the solutions in the system

{P (x) = 0 ∈ R n-1 , x 1 -α = x 2 -β = 0} is at most two, i.e.,
B := pop(L) 4:
if 0 ∈ P (B) and rank( J P (B)) < n -1 then 5:

Subdivide B and add its children to L.

6: return True.

Lemma 5.2.3. Semi-algorithm 1 stops if and only if P satisfies Assumption A 1 in B.

Proof. If Semi-algorithm 1 stops, by the conditions in Step (4), the box B is partitioned into two sets of boxes. A set of boxes that are disjoint with C and the other one is a set of boxes that contain parts of C that satisfy Assumption A 1 . Thus, Assumption A 1 is satisfied in B. On the other hand, assume that P satisfies Assumption A 1 in B and Semi-algorithm 1 does not stop, then, for every positive real there exists a closed box B ⊂ B, with w(B ) < such that the conditions in Step (4) are satisfied in B . Consider the infinite chain

B 1 1 , B 1 2 , B1 3
. . . and take q k ∈ B 1 k , with q k = q k for k = k . Since B is compact, then there exists a subsequence of q k that converges to a point on B say q. Since P and J P are box function we must have that P (q) = 0 and rank(J P (q)) < n -1. Thus, q is a point in C that does not satisfy Assumption A 1 which is a contradiction. Hence, Semi-algorithm 1 stops. The next step is to check Assumptions ℵ 1 and ℵ 2 . For this goal, we want to find a finite set of pairwise disjoint boxes in B Ball such that every box contains at most one solution of Ball(P ) and the union of these boxes contains all solutions of Ball(P ) in B Ball . Notice that, by the definition of box functions, for a closed (2n-1)-box U, if 0 ∈ Ball(P )(U), then U does not contain a solution of Ball(P ), whereas the condition 0 ∈ Ball(P )(U) does not necessarily imply that a solution is in U. This is why the set we are going to find might have unnecessary boxes. However, we will see later that this is enough for our purpose. Before introducing Semi-algorithm 2, we define the following functions.

Definition 5.2.4. Consider the set R 2n-1 t 0 = {(x 1 , x 2 , y, r, t) ∈ R × R × R n-2 × R n-2 × R | t 0} and define f ± Ball : R 2n-1 t 0 → R n (x 1 , x 2 , y, r, t) → (x 1 , x 2 , y ± r √ t) Define the function f Ball : R 2n-1 t 0 → R n × R n that maps X to (f + Ball (X), f - Ball (X)).
Notice that f Ball is an extension of Ω P (Definition 4.1.7). By abuse of notation, for a set S ⊂ R 2n-1 , we define f Ball (S) as f Ball (S ∩ R 2n-1 t 0 ).

Semi-algorithm 2: Isolating the solutions of Ball(P ) (under Assumption A 1 )

Input: An open n-box B, a function P from B to R n-1 such that P satisfies Assumption A 1 in B and a (2n -1)-open box U 0 that contains B Ball (see Remark 5.2.7). Termination: If and only if Ball(P ) satisfies Assumptions ℵ 1 and ℵ 2 in B Ball . Output: A list of pairwise disjoint isolating boxes of the solutions of Ball(P ) in U 0 such that their images by f Ball lies in B × B.

1: Solutions = ∅. 2: L := {U 0 }. 3: while L = ∅ do 4: U := pop(L). 5: if 0 ∈ Ball(P )(U) or ( f Ball (U)) ∩ (B × B) = ∅ then 6:
Do nothing (U is simply removed from L).

7:

else if rank( J Ball(P ) (U)) = 2n -1 and f Ball ( -inflation(U)) 6 ⊂ B × B then 8:
if -inflation(U) contains a solution of Ball(P ) (see Remark 5.2.5) then 9:

Add -inflation(U) to Solutions. Subdivide U and add its children to L.

12: Remove duplicates in Solutions (see Remark 5.2.5). 13: return Solutions Remark 5.2.5. Steps (8) and ( 12) are not detailed because they are standard in subdivision algorithms to handle the issue of solutions on or near box boundaries and ensuring that solution boxes are pairwise disjoint. We only sketch below how these steps are done and refer to [Sta1995, §5.9.1 ; Kea1997 ; XY2019] for details. In Step (8), an existence test is performed by evaluating an interval Newton operator on an -inflation of the box U. The inflated box -inflation(U) is certified to contain a solution if its image by the interval Newton operator is contained in the interior of -inflation(U). When the existence test is positive, the solution may be on the boundary or even outside U, but still in the interior of -inflation(U). The side effect is that the same solution may be reported several times when it is on or near a boundary of the subdivision. This issue is then solved in Step (12) as follows. When two boxes in the set Solutions intersect, they must report the same solution, and in addition, this solution is in the intersection of the two boxes. In

Step (12), we thus compute intersections between boxes and replace intersecting ones by their intersection box. The boxes in the output set is Solutions are thus pairwise disjoint.

Lemma 5.2.6. Under Assumption A 1 in B, if Semi-algorithm 2 stops, it returns a list of pairwise disjoint isolating boxes of the solutions of Ball(P ) in U 0 such that their images by f Ball lies in B × B. Moreover, Semi-algorithm 2 stops if and only if Ball(P ) satisfies Assumptions ℵ 1 , ℵ 2 in B Ball .

Proof. We first prove the correctness of the Semi-algorithm 2 assuming that it terminates. Since Step (5) is the only time the algorithm discards boxes, it never discards a box that contains a solution of Ball(P ) in U 0 such that its image by f Ball lies in B × B. Hence, all such solutions of Ball(P ) lie in output boxes. The rank condition in Step (7) guarantees that each output box contains at most one solution of Ball(P ) [Sny1992, Theorem A.1]. The fact that every output box contains at least one solution is ensured by a standard algorithm in Step (8) (see e.g., [Neu1991, Theorem 5.6.2 ; XY2019] and Remark 5.2.5). Finally, by Step (12), the output boxes are pairwise disjoint, hence the algorithm outputs isolating boxes of the solutions of Ball(P ) in U 0 such that their images by f Ball lie in B × B.

To prove the equivalence for the termination, first assume that Semi-algorithm 2 stops and returns Solutions. According to the correctness proof, every solution X of Ball(P ) in B Ball is regular and satisfies Ω P (X) ∈ B × B. Thus, Assumptions ℵ 1 and ℵ 2 are satisfied in B Ball .

On the other hand, assume that ℵ 1 and ℵ 2 hold in B Ball . We prove that Semialgorithm 2 terminates. By Assumption ℵ 1 all solutions in B Ball of the square system Ball(P ) are regular. Hence, they form a zero dimensional manifold in the compact space B Ball . Thus, the solution set is finite. We now prove that for any box U ∈ L with a small enough width, one of the conditions in Step (5) or the conditions in Steps (7-8) are satisfied. Thus, in both cases U will be removed from L, and hence, Semi-algorithm 2 stops after a finite number of iterations. Due to Assumption ℵ 2 , after a finite number of iterations, no box U in L intersects the boundary of B × B. Moreover, due to the convergence of the box evaluations, we can also assume that either f Ball ( -inflation(U)) ⊂ B × B, which is the second condition of Step (7), or ( f Ball (U)) ∩ (B × B) = ∅ which is the second condition of Step (5).

If U does not contain a solution of Ball(P ), then due to convergence of the box function evaluation of Ball(P ), after a finite number of iterations, every children

U of U satisfies 0 ∈ Ball(P )(U ), that is, it is discarded in Step (5).
If U contains a solution of Ball(P ) in B Ball , according to Assumption ℵ 1 , it is a regular solution. Due to the convergence of the box evaluation det(J Ball(P ) (U)) will eventually be non zero and thus rank( J Ball(P ) (U)) will eventually be 2n -1 after a finite number of iterations, which is the first condition of Step (7). Due to the convergence of the interval Newton existence test, the condition of Step (8) will also be eventually satisfied (see Remark 5.2.5). The refined box will then eventually be added in the Solutions list.

Thus, for any box in L with a small enough width, one of the conditions of Step (5) is satisfied or all of the conditions in Step (7-8) are satisfied, thus it is either discarded or added to the output. Hence, Semi-algorithm 2 terminates.

Remark 5.2.7. Semi-algorithm 2 requires a closed (2n -1)-box U 0 that contains B Ball . For instance the following set could be used:

{(q, r, t) ∈ R 2n-1 | q ∈ B, -1 r i 1 for i ∈ {3, . . . , n}, 0 t ξ 2 4 } with ξ = max { q -q | q, q ∈ B}.
Finally, using Lemma 5.2.2, Semi-algorithm 3 checks whether P satisfies Assumptions A 1 , A 2 , A 3 , A 4 and A - 5 in B and outputs a superset of isolating boxes of the singularities of π C (C). Output: N , a list of certified node singularities: a list of boxes in R 2n-1 whose projections in R 2 contain each a single node of π C (C). U , a list of uncertified singularities: a list of boxes in R 2n-1 whose projections in R 2 contain each at most one node or one cusp of π C (C).

The union of all these projected 2D boxes contains all the singularities of π C (C).

1: Check Assumption A 1 (Semi-algorithm 1). 2: Compute a closed (2n -1)-box U 0 that contains B Ball (Remark 5.2.7).

3: L := output of Semi-algorithm 2. 4: Keep refining all boxes U ∈ L (see Remark 5.2.8) until no triplets of boxes overlap in projection. Then remove from L one box from every pair (see Remark 4.1.10). This ensures Assumption ℵ 3 .

5: N := boxes of L that lie in the halfspace t > 0.

6: U := boxes of L that intersect the hyperplane t = 0. 7: return N and U .

Remark 5.2.8. The refinement of an isolating box of a solution is performed by iterative evaluation of an interval Newton operator; we refer to [Neu1991, Theorem 5.6.2] for details.

To identify the possible cusp points in the set U returned by Semi-algorithm 3, one may wish to solve independently the Ball system with the additional constraint t = 0 (by Remark 4.1.9). Unfortunately, in this case we have an over-determined system and thus we cannot certify its solutions numerically. In the special case of a silhouette curve, it is possible to identify cusp points with numerical algorithms in the case n = 3 [IMP2016a, Lemmas 9 & 10], but we leave as a conjecture its generalization for n > 3.

On the other hand, for curves that satisfy the strong assumptions, A 5 ensures that there are no cusps in the projection, which is equivalent to L c being empty and Ball(P ) having no solutions on the hyperplane t = 0 (by Remark 4.1.9). Hence, if Assumptions A 1-5 hold, we can refine the boxes output by Semi-Algorithm 3 until no box intersects t = 0. Boxes in the half-space t < 0 correspond to imaginary points in C n (Definition 4.1.7). Then the boxes satisfying t > 0 are all the isolating boxes of the nodes of π C (C) by Lemmas 5.2.6 and 4.1.6, Remark 4.1.9 and Corollary 4.3.5. Termination: If and only if P satisfies Assumptions A 1-5 in B. Output: A list of boxes in R 2n-1 whose projections in R 2 are isolating boxes of the singularities π C (C) (all singularities are in some boxes and each box contains a unique singularity).

Semi-Algorithm 4: Checking Assumptions

1: N, U := output of Semi-Algorithm 3.

2: for all U ∈ U do 3:

Keep refining U (see Remark 5.2.8) until it does not intersect the hyperplane t = 0 and discard it if it lies in the half-space t < 0.

4: return N ∪ U .

Semi-algorithms improvements

Semi-algorithm 3 (resp. 4) of Section 5.2 takes as input an open n-box B and a curve C defined by a smooth function P from B to R n-1 . It terminates if and only if P satisfies our weak (resp. strong) assumptions of Definition 3.1.1. Upon termination, it outputs a superset of the singularities (resp. the singularities) of the curve π C (C). This section is dedicated to describe the (semi) algorithm we implemented which is a variant of Semi-algorithms 3 and 4.

In the following, we present the challenges we faced and how we overcame them for implementing Semi-algorithms 3 and 4.

Subdivision-based solver on a high-dimensional system

The first challenge we faced is solving a high-dimensional system such as the Ball system which is of dimension 2n -1. The solutions of the Ball system are detected in Semi-algorithm 3 using a method based on the subdivision of the space. Such methods usually suffer from the curse of dimensionality in practice. In the following, we describe the way we overcame this challenge.

For the rest of Section 5.3.1, we assume that Assumption A 1 is always satisfied. Moreover, recall that M i is the sub-matrix of J P obtained by removing the i-th column. Over a closed n-box B, we call P parameterizable in x i if M i is full rank at every point of B. This notion has a geometric interpretation: the function P being parameterizable in x i implies that the tangent space at any point, of C, is not orthogonal to the x i axis.

The domain B Ball where the Ball system is solved is refined by first enclosing the curve C in the smaller space R n . Our approach follows the observation that if the curve is enclosed by a set enclosing_curve of boxes in R n such that:7 (i) The intersection of the plane projections of two boxes in enclosing_curve. These 2-boxes contain the nodes induced by the projections of two points in C, with the same projection, that do not belong to one and the same box of enclosing_curve.

(ii) The projection of a box B in enclosing_curve such that the curve C in this box is not parameterizable in x 1 and not parameterizable in x 2 . Such a box can be detected using the condition 0

∈ det( M 1 (B)) ∩ det( M 2 (B)),
where M 1 and M 2 are box functions of M 1 and M 2 respectively. These 2-boxes contain the nodes induced by the projections of two points, with the same projection, that belong to one and the same box B of enclosing_curve. Figure 5.1 illustrates an example of this case where B contains a part of C that projects to a loop. In this figure, to give a geometric intuition, we presume that q 1 and q 2 are in C.

Moreover, every cusp of π C (C) lies in the projection of a box of enclosing_curve containing a point of C with a tangent orthogonal to the (x 1 , x 2 )-plane, that is, P is parameterizable neither in x 1 nor in x 2 at the same point in C. These boxes are contained in the set of boxes that satisfy the constraint (ii) above.

The above observation helps to cover the set L c ∪ L n using enclosing_curve. This paves the way to cover all Ball solutions using Lemma 4.1.6. We are going to consider this approach to prune B Ball which makes Semi-algorithms 3 and 4 more efficient in practice. However, the reader might notice that any two neighbor boxes in R n satisfy the condition in (i). Hence, a non-ignorable set of spurious pairs satisfy (i). This makes the previous condition a bad criteria to detect points of L n . In the following lemma, we develop the previous observation to obtain the criteria by which we cover L c ∪ L n with less redundancy than the above observation.

Lemma 5.3.1. Let C be a curve defined by P satisfying Assumption A 1 in a closed n-box B. Then, every point of L c ∪ L n is in at least one of the following: mance of the algorithm. Item (c) is to preserve the connectivity of the components of C. This will be used in Lemma 5.3.1 where the connectivity of C in a box plays a role in detecting points of Ln in our approach.

x 1 x 2 (x 3 , . . . , x n )

q 1 q 2 B ⊂ R n C π C (C) Node Figure 5.1: An example of Case (ii) for a curve C in R n .
The box B contains a smooth part of C that projects to a loop with a node. This implies that two points q 1 and q 2 are in B such that det(M 1 (q 1 )) = 0 and det(M 2 (q 2 )) = 0. Hence, the tangent space plane projection of q 1 (resp. q 2 ) is parallel to x 2 -axis (resp. x 1 -axis).

(1) The union of disjoint boxes B and B in enclosing_curve such that their plane projections intersect.

(2) The union of two non-disjoint boxes B and B in enclosing_curve such that C has at least one x 1 -critical point and at least one x 2 -critical point in B ∪ B . Proof. For a points q ∈ L c , since T q C is orthogonal on R 2 , then q is x 1 -critical and x 2 -critical. Hence any box B in enclosing_curve, containing q, satisfies (2) (taking B = B). On the other hand, let q in L n and q be a point in L n with π C (q) = π C (q ). Suppose that q and q are in (possibly equal) boxes B and B respectively. We distinguish the following cases:

(a) B ∩ B = ∅: Then, Case (1) applies; recall that the projections of B and B intersect since both contain π C (q) = π C (q ). q π(q ) = π(q ) q q π(q B and B at q and q respectively. The curve and the boxes intersect at q such that q does not have the same plane projection as q and q . We are going to prove that B and B satisfy (2) using Lagrange multiplier theorem (LMT). We consider x 1 as a function from R n to R and P to be the constrains function. Since C is connected in B ∪ B and since π C (q) = π C (q ), there must be an x 1 -critical point q in the compact subset C in B ∪ B .8 By LMT, there exists a non-trivial λ ∈ R n-1 such that J x 1 (q) = λ t J P (q). This implies that λ t M 1 (q)) = 0 ∈ R n-1 . Thus, det(M 1 (q)) = 0. In other words, q is a x 1 -critical point. In the same way, we prove the existence of q in C restricted to B ∪ B such that det(M 2 (q )) = 0. This proves that q is a x 2 -critical point. Thus, B and B satisfy (2).

) B B C π C (C)
Remark 5.3.2. From the proof of Lemma 5.3.1, we can obtain a more precise statement: Every point of L c is in a box B of enclosing_curve that satisfies (2), considering B = B. Let q and q be in C with π C (q) = π C (q ), in particular, q and q are in L n . Assume that q and q are in the (possibly equal) boxes B and B respectively. If B and B are disjoint, then they satisfy (1). Otherwise, they satisfy (2).

An interesting question is how to detect boxes satisfying (1) or (2). Detecting boxes of enclosing_curve satisfying (1) can be done using any black-box tool that reports the pairs of intersected boxes of a given set of boxes. See Remark 6.1.1 for information about the tool we chose. Boxes B, B satisfying (2) can be detected using the same tool in addition to checking the condition:

0 ∈ M 1 (B) ∪ M 1 (B ) and 0 ∈ M 2 (B) ∪ M 2 (B ).
The next step is to cover the Ball solutions. This is be done by the following proposition which is a corollary of Lemma 5.3.1.

Proposition 5.3.3. Let C be covered by enclosing_curve as above. Let X ∈ B Ball be a solution Ball(P ). Then, at least one of the following holds:

(a) There exists B in enclosing_curve which is the cross product of boxes Proof. Let X = (x 1 , x 2 , y, r, t) ∈ R×R×R n-2 ×R n-2 ×R be a solution of Ball(P ) in B Ball . Let (q, q ) be a pair of L (Definition 4.1.5) such that f Ball (X) = (q, q ). Recall that q and q are in L c ∪L n . Hence, by Remark 5.3.2, q and q are respectively in two boxes B and B of enclosing_curve satisfying (1) or (2). In both cases, we imply that X ∈ f -1

B (x 1 ,x 2 ) in R 2 and B y in R n-2 satisfying (2) and X lives in f -1 Ball (B × B). Hence, X is in (x 1 , x 2 , y, r, t) ∈ R × R × R n-2 × R n-2 × R, where (x 1 , x 2 , y) ∈ B, r ∈ [-1, 1] n-
(2) such that X is in f -1 Ball (B × B ). Hence, X lives in (x 1 , x 2 , y, r, t), with (x 1 , x 2 ) = B (x 1 ,x 2 ) ∩ B (x 1 ,x 2 ) , y = f -1 Ball B Ball R 2n-1 C B ⊂ R n
Ball (B × B ) which proves the proposition.

We now reformulate the results of Lemma 5.3.1 and Proposition 5.3.3 in Algorithm 5. This algorithm is used to find a set of (2n -1)-boxes that cover the Ball solutions instead of solving it in the whole B Ball . It turns out in practice, that computing this set and solving Ball(P ) over its boxes is faster than solving the Ball system over B Ball .

Algorithm 5: Pruning B Ball .

Input: An open n-box B and a curve C defined by P ∈ C ∞ (R n , R n-1 ), with C ⊂ B such that Assumption A 1 is satisfied. Output: A finite set Solving_domains which is a set of (2n -1)-boxes contained in B Ball that contains all the solutions of Ball(P ) in B Ball . )]

6:

Add (x 1 , x 2 , y, r, t) to Solving_domains. Add (x 1 , x 2 , y, r, t) to Solving_domains.

7: for distinct B = B (x 1 ,x 2 ) × B y and B = B (x 1 ,x 2 ) × B y in enclosing_curve satisfying (1) or (2) do 8: (x 1 , x 2 ) = B (x 1 ,x 2 ) ∩ B (x 1 ,x 2 ) and y = 1 2 (B y + B y ) 9: r = [-1, 1] n-2 if B ∩ B = ∅,

11: return Solving_domains

As a part of the improvements, the output of Algorithm 5 replaces L in Semi-Algorithm 2, Step (2). For more details, Section 5.3.3 is dedicated to provide a global vision of the improvements.

Evaluating the operator D

To solve the Ball system we need box functions for the operators S and D (Definition 4.1.1). We first note that if P (x 1 , x 2 , y) ∈ C ∞ (R n , R) is a polynomial function, S • P (resp. D • P ) is also a polynomial function and can be computed from the terms of P (x 1 , x 2 , y + rt) that have even (resp. odd) exponents in the variable t, see [IMP2016b, Lemma 6] for details.

If P is a more general C ∞ function for which we have a box function, using Definition 4.1.1, computing S • P on any (2n -1)-box U or computing D • P on U such that its t-interval does not contain 0, is implemented from the box function of P and interval arithmetic. On the other hand, when the t-interval contains 0 or is close to 0, the division by √ t in the formula for D • P makes the computation undefined or unstable. In such non-polynomial cases, we use a Taylor form at order 3 [Ral1983], that is, we compute a Taylor expansion with remainder at t = 0 of D • P and evaluate by interval the third order derivative. We define a threshold δ Taylor (which we set to 10 -2 ) such that this Taylor form is used when the t-interval has values smaller than δ Taylor .

Improvement integration

In this section, we discuss how we modify Semi-algorithms 2, 3 and 4 in view of the improvements we discussed in Sections 5.3.1 and 5.3.2. We also add a new parameter eps min in addition to changes in several steps in Semi-algorithms 2, 3 and 4, preparing the semi-algorithms to the implementation stage.

Parameter eps min . We start with adding the parameter eps min to Semi-algorithm 4 and any (semi) algorithm called by it. The parameter eps min is a non-negative real that guarantees the termination in Semi-algorithm 4 even if the termination criteria is not satisfied. More precisely, if eps min = 0, then the parameter does not affect Semi-algorithm 4. Otherwise, this parameter forces the termination of subdivision or refinement steps in Semi-algorithm 4, that is, the semi-algorithm does not subdivide (or refine) further if the width of the processed box (which is a child of B) is less than eps min . In this case, the remaining boxes (where eps min forces the termination of subdivision or refinement) are labeled as "unknown-status". These boxes are left aside and not considered for the further steps. For example, in Semialgorithm 2 (which is called by Semi-algorithm 4), Step (11), the parameter eps min prohibits the subdivision when the width of U is less than eps min . From a practical point of view, the parameter eps min allows the use of fixed precision arithmetic. 9In addition, it is useful if we want to ensure the termination even when not all assumptions A 1-5 are satisfied, certifying π C (C) singularities only on the part for which we verify the strong assumptions. On the other hand, an "unknown-status" label of a box prohibits the use of the Ball system to certify the singularities of its plane projection (if any). This is, of course, a drawback of setting eps min > 0.

Changes in Semi-algorithm 4. We now describe the changes in Semi-algorithm 4 and all (semi) algorithms called by it. i. Assumption A 1 is not checked by Semi-algorithm 1 anymore. It is rather checked using a black-box solver (see Black-box solver 7 in Section 6.1), considering the parameter eps min mentioned above. Such a solver does not only check Assumption A 1 , but also computes a set enclosing_curve of nboxes that enclose C with Properties (a), (b) and (c) in Section 5.3.1. More precisely, for eps min > 0, the solver returns possibly-empty two sets of nboxes. The first enclosing_curve encloses a part of C such that the solver succeeds to verify Assumption A 1 . The second set contains boxes of width less than eps min and encloses the rest of C.

ii. The set enclosing_curve is required for Algorithm 5. Furthermore, the output of Algorithm 5 is used instead of U 0 in Step (2) of Semi-algorithm 2. Consequently, Step (2) of Semi-algorithm 3 is omitted. In practice, solving the Ball system over Algorithm 5 output is more efficient than solving over U 0 , despite the fact that computing U 0 is simpler.

iii Considering the parameter eps min and the changes above, we reformulate the specifications of Semi-algorithm 4:

Semi-Algorithm 6: Variant of Semi-algorithm 4

Input: An open n-box B, a smooth function P from B to R n-1 and two nonnegative reals δ Taylor and eps min .

A 1-5 in B or eps min > 0. Output: N , a list of pairwise disjoint (2n -1)-boxes whose projections in R 2 are isolating boxes of node singularities in π C (C). U , a list of n-boxes and (2n-1)-boxes such that Assumption A 1 is not verified (n-boxes), or Assumptions ℵ 1 , ℵ 2 or ℵ 3 are not satisfied ((2n -1)-boxes). The width of every box is less than eps min . Moreover, every singularity is in the plane projection of a box of N or U .

In the next chapter, we provide more details about the implementation and how to use the software we develop based on Semi-algorithm 6 and the improvements explained in this section. This part details our implementations of the (semi) algorithms presented in Chapter 5 and reports on our experimentations. Based on the theoretical arguments in Chapters 3 and 4 and the (semi) algorithms presented in Chapter 5, we developed the software Isolating_singularities11 written in Python. Given a curve C in R n , our software finds the singularity of its plane projection. In addition, the software provides a tool to visualize the projected curve and its singularities. Chapter 6 is dedicated to describe Isolating_singularities in more details.

In chapter 7, we examine the efficiency of our software performing four experiments that show different cases: C is an analytic curve in R 3 or R 4 (Sections 7.1 and 7.2), C is a high-degree algebraic curve defined by randomly-generated polynomials (Section 7.3) and C is a curve that contains two very close branches whose projections define a node (Section 7.4).

Finally, in Chapter 8, we revisit a classical problem in robotics design and analysis. This problem consists of determining the configurations in which a robot is allowed to move without damaging it. More formally, it consists of numerically computing the generalized aspects of that robot (Definition 8.1.1 (g)). For a robot of two degrees of freedom, the boundaries of the generalized aspects is generically a space curve in R 4 . We use our software on the fully parallel robot RRRRR and we present our experiments in Section 8.3. We present in this chapter a brief description of our software Isolating_singularities12 which is an implementation of Semi-algorithm 6.

Our software consists of three main functions:

enclosing_curve, enclosing_singularities and ploting_boxes. The first one encloses the curve C in R n and checks Assumption A 1 . More precisely, it returns (a) a set of n-dimensional boxes that cover a smooth part of the curve C, and (b) a set of n-dimensional boxes that cover an unknown part of C, where the solver fails to certify the smoothness.

The goal of the second function enclosing_singularities is to isolate the Ball system solutions. Consequentially, the plane projections of the boxes output by enclosing_singularities are, in the best case,13 isolating boxes of the singularities of the projected curve, and in the worst case, they enclose the singularities. More precisely, the projection of the boxes returned by enclosing_singularities consists of:

• a set of 2-dimensional boxes each of which contains a unique singular point of π C (C), and

• a set of 2-dimensional boxes where the behavior of π C (C) is unknown, but such that every singular point in π C (C) is contained in a box of one of these two sets.

The third function ploting_boxes performs a visualization task. Using the Python Package Matplotlib [Hea2007], we provide a tool to 2D-visualize the results of the first and second functions.

In Section 6.1, we present the third-party libraries that we used in our software. The second part of this chapter (Section 6.2) documents the software Isolating_singularities.

Third-party libraries

Our software is based on interval arithmetic, interval evaluations of analytic functions and an interval solver. We use the following libraries, Python-FLINT [Joh2012b] and Ibexsolve [Nin2015] for these tasks.

Python-FLINT is a Python extension module wrapping FLINT (Fast Library for Number Theory) and Arb (arbitrary-precision ball arithmetic) [Joh2012a], which offers a toolbox for interval arithmetic and evaluation of analytic functions.

Ibexsolve is a C ++ end-user program that solves systems of non-linear equations rigorously, that is, it does not lose any solution and returns each solution under the form of a small box enclosing the true value. It implements a classical branch-andprune algorithm that interleaves contractions and branching (bisections) to enclose the solutions of a system at any given desired precision. However, as opposed to Arb, Ibexsolve has a fixed precision, hence when several solutions are closer to each other than this precision, it will correctly return an enclosing box for these solutions but it will fail at isolating them. In our software, we use the default precision which is 10 -7 . Ibexsolve, and thus also our software Isolating_singularities, use a parameter eps_max that defines a maximum width for the isolating boxes output by the solver (the box bisections are forced until all output boxes are not larger than eps_max). We use Ibexsolve for solving the Ball system (Semi-algorithm 2) and also in a variant of Semi-algorithm 1 to check the smoothness of the curve C and at the same time enclosing C in a set of boxes of R n . The following is a specification of Ibexsolve:

Black-box solver 7: Ibexsolve Input: A system of equations, the ranges of the variables and non-negative reals eps min eps max . 14 Output: The following sets of boxes whose union contains all solutions of the given system: (a) A set of certified boxes, of width at most eps min , where every box contains a regular solution (resp. full-ranked-Jacobian part of the solution set) if the system is zero-dimensional (resp. of positive dimension).

(b) A set of boxes, of width at most eps max , that intersect the boundary of the input box where the variables are defined and contain solutions of the given system. Ibexsolve provides an option to verify the regularity of the solutions inside these boxes.

(c) A set of unknown-status boxes of width at least eps min .

Remark

The Python software Isolating_singularities

The following is a short specification of the software Isolating_singularities16 which is a Python implementation of Semi-algorithm 6. As we mentioned in the introduction of this chapter, two main tasks can be performed using our software: computing the output of Semi-algorithm 6 (the function enclosing_singularities) and a visualization tool (ploting_boxes). There is an additional function enclosing_curve that performs a sub-task, namely, computing a set of n-boxes that enclose the curve C. The output of this function is a part of the previous functions input. The following is a more detailed description of how to use these functions. By an interval we mean a list [a, b] of floats a, b that represent the lower and upper bounds respectively. By a box we mean a list of intervals.

• enclosing_curve(system, Box, X, eps_min, eps_max):

This function returns a list of two lists: The first one is a list of boxes in R n that cover a smooth part of the curve (more precisely, where the used solver Ibex-Solve succeeds to certify the smoothness). The second one contains boxes where IbexSolve fails to certify the smoothness of the curve.

Parameters.

system: a string that is the name of the txt file which contains the equations that define the curve

-Box: a box in R n
-X: a list of the SymPy symbols that appear in the equations of system -eps_max : the maximal width of the output boxes, that is, this is a criterion to force bisection: the width of output box will not be larger than eps_max. By default, eps_max = 0.1.

-eps_min : the minimal width of output boxes. This is a criterion to stop bisection: for a box where IbexSolve fails to certify the regularity of the solutions inside that box, the width of this box will not be larger than eps_min. By default, eps_min = 10 -7 .

• enclosing_singularities (system, enclosing_curve, Box, X, eps_min, eps_max, delta_Taylor )

The output of this function consists of the following two lists such that every solution of Ball(P ) (no matter it is regular or not) lives in the union of these sets of boxes:

(a) A set of enclosing boxes, of width at most eps max , each of which contains one certified (regular) solution of Ball(P ). Thus, the plane projection of every box contains a unique singularity, namely, a node.

(b) A set of boxes of widths at most eps min such that, in their projections, the behavior of the plane curve is unknown.

Parameters.

system, Box, X, eps_min and eps_max: are as in enclosing_curve.

-enclosing_curve: is the output of function enclosing_curve.

threshold: Represents the parameter δ Taylor introduced in Section 5.3.2. When evaluating the Ball system, the parameter controls when to use the system in (4.1) or the Taylor form. The default value of this parameter is 10 -2 .

• ploting_boxes(certif ied_boxes, uncer_boxes, var, Box, nodes)

Based on NumPy [HMvdW + 2020] and Matplolib [Hea2007], this function plots a plane projection of enclosing_curve output and enclosing_singularities output.

Parameters.

certif ied_boxes : A list of boxes that cover a smooth part of the curve certified by IbexSolve. The function plots the plane projection of this list in green.

-uncer_boxes : A list of boxes that cover parts of the space curve with unknown status. The function plots the plane projection of this list in red.

var (optional): A list of two integers that determines the variables for which the plane projection is considered. By default, the function considers the first two variables.

- In this section, we present four experiments performed with our software Isolating_singularities. The first example is pedagogical and considers a simple analytic curve in R 3 that induces only one node and one ordinary cusp in R 2 . The second example considers a smooth analytic curve in R 4 that induces many nodes in R 2 . The third one considers sparse but reasonably-high-degree algebraic equations in R 4 . It should be stressed that, up to our knowledge, the two latter examples are out of reach by other methods: indeed, no other certified algorithm can handle non-algebraic curves in dimension higher than 3 and, for reasonably-high-degree algebraic equations in R 4 , the bivariate equation defining their 2D projection often has a very high degree (see Section 7.3 for details).

Finally, in the fourth example, we exhibit the behavior of our software when a node in R 2 is induced by a pair of points (on the space curve) that are very close. Indeed, when the equations defining the space curve are not algebraic, the Ball system contains a division by √ t (due to the formula of D • P ), which may cause instability since t tends to zero when the distance between the pair of points tend to zero. For that purpose, we consider two very close skew lines defined analytically (and not algebraically).

We report the running times and other relevant parameters in Tables 7.1 and 7.2. Running times are in seconds on a sequential Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz machine with Linux. We emphasize that the experiments are done with a prototype implementation that is under ongoing development. The tree size columns reports the total number of boxes created during the subdivision algorithm either for enclosing the curve C in R n or for solving the Ball system in R 2n-1 . For the enclosing part, the column output boxes is the number of boxes in the set enclosing_curve. For each experiment, we provide a visualization of the plane projected curve π C (C) with its singularities. On each figure, the green boxes are the plane projections of the boxes in enclosing_curve that enclose C, hence these green boxes enclose π C (C). The black boxes are the projections of the Ball system solution boxes identifying nodes of the plane curve π C (C).

For each experiment, we consider three values of eps_max and it can be observed (see Table 7.1) that the smaller the value of eps_max, the larger the set enclosing_curve, and the longer it is to compute. As expected, even with the improvement to reduce the Ball system domains to be solved in, the subdivision in the high dimensional space R 2n-1 is the dominant step of the algorithm.

Experiment 1: Analytic curve in R 3 generating one node and one ordinary cusp

We start with a pedagogical example pictured in Figure 7.1. Running times are given in Table 7.1. The curve C is defined in the box B = (-1, 4) × (-1, 4) × (-4.8, -1.4) by

P (x 1 , x 2 , x 3 ) = [x 1 -cos(x 3 )(3 + sin 4 (x 3 )) + 3, x 2 -sin 2 (x 3 )(3 + sin(2x 3 ))].

Plane projection of boxes in enclosing curve

Plane projection of Box2 Plane projection of Box1

x 1

x 2

Figure 7.1: Experiment 1: Plane projection of an analytic curve in R 3 with one node and one ordinary cusp.

Our improved Semi-algorithm 3 outputs the following solutions for the Ball system in R 5 : Box 1 in the set N thus projects to a node of π C (C). Box 2 being in the set U , one cannot decide whether its projection in the plane contains a node, a cusp or no singularity at all. On the other hand, one can notice on the equation P = 0 that the curve is parametrizable by the variable x 3 . It is thus an easy computation to check that for the value x 3 = -π, the point q = (0, 0, -π) is on the curve C and its tangent line at q is generated by the vector (0, 0, 1) which is orthogonal to the projection plane. It is then clear that the projection of q generates a cusp that is witnessed by Box 2 .

N = {Box1 = [3, 3] × [3, 3] × [-3.15, -3.14] × [1, 1] × [2.
7.2 Experiment 2: Analytic curve in R 4 with many nodes

P = [x 1 + 2 sin(x 1 ) -cos(x 4 ) -(3 cos(x 3 ) -cos(2.8571x 3 )),
x 2 + 0.2 cos(x 2 ) + (3 sin(x 3 ) -sin(2.8571x 3 )) + sin(x 4 ),

x 2 4 -sin(x 2 )] over the box B = (-1, 0) × (-0.1, 3.5) × (-20, 20) × (-10, 10). This curve has many nodes, some of them very close to each other. Running times are given in Table 7.1.

Experiment 3: High degree algebraic curve in R 4

The goal of this experiment is to emphasize the genericity of the assumptions and the efficiency of our software in the sparse polynomial case. The curve C in defined in the 4-box B = (-1, 0.2) × (-0.2, 1.4) × (-10, 10) 2 and is the zero set of three polynomials of degrees 17, 15 and 13, respectively that have a unique monomial of highest degree (which is monic) and 9 other random monomials of degrees at most 2 with integer coefficients in (-25, 25). using elimination theory and its singularities can then be isolated using algebraic solvers. However, the implicit equation we obtained for π C (C) is defined by an irreducible bivariate polynomial of degree 442 with 51074 monomials. Isolating the singularities of such a high-degree polynomial is then a real challenge.

P = [x 17 1 -14x 2 1 -7x 1 x 3 -7x 2 2 -22x 2 x 4 -x 3 x 4 -19x 2 4 + 8x 1 -14x 3 + 9, x 15 2 + 8x 1 x 3 -14x 1 x 4 -15x 2 2 + 16x 2 x 3 + 8x 2 x 4 + 2x 2 3 + 13x 2 4 + 11x 1 + 11x 2 , x 13 3 + 17x 2 1 -15x 1 x 2 + 4x 1 x 3 -20x 1 x 4 + 2x 2 2 -10x 2 x 3 + 4x 2 x 4 + 20x 2 4 -23x 2 ].
Note also that our class of examples is rather specific and our software does not work that well if our defining polynomials are dense or even if they have nonmonic high-degree monomials. However, it should be stressed that our software is a prototype and that it would be more efficient to use an interval solver specialized for the algebraic case than the versatile Ibex solver we used. Such a specialized approach has been proved successfully for generic curves in the 3-dimensional case by Imbach et al. [IMP2018].

Experiment 4: Two close lines in R 3 generating a node

As mentioned in the preamble of Chapter 7, the purpose of this experiment is to study the behavior of our software when a node in R 2 is induced by a pair of points (on the space curve) that are very close; namely when a node (x 1 , x 2 ) ∈ R 2 is induced by the pair of points (x 1 , x 2 , y ± r √ t) ∈ R n with t that tends to zero. Indeed, when the equations defining the space curve are not algebraic, the Ball system contains a division by √ t (due to the formula of D • P ), which may cause instability since t tends to zero when the distance between the pair of points tend to zero.

The simplest example to consider is the two skew lines x 2 = x 1 in the plane x 3 = and x 2 = -x 1 in the plane x 3 = -, whose projection in the (x 1 , x 2 )-plane has a node at the origin, and to make vary towards 0. The pair of lines is thus defined by [ x 2 -x 1 x 3 , (x 3 -)(x 3 + )] but, in order to have non-algebraic equations, we replace x 3 by sin x 3 and define our two lines by

P = [ x 2 -x 1 sin x 3 , sin 2 x 3 -2 ] in the box B = (-1, 1) 3 .
The goal of this experiment is to illustrate the stability of our software when varies towards 0. Recall from Section 5.3 that the D operator is evaluated on a box in two different ways depending on how close to zero is the t-interval of that box. The Ball system is thus solved either with Equation (4.1) (involving a division by √ t) when the values of the t-interval are larger than a parameter δ Taylor (set to 10 -2 ) or using a Taylor expansion otherwise.

To illustrate the stability of our software, we compared in Table 7.2 its running times when varies towards 0 with what it would be without using Taylor expansions. It shows that if we do not use Taylor expansions, the solution is not certified (by the Ibexsolver; see Section 6.1) when 10 -5 . On the other hand, our software is stable, although its running time increases from 0.1 to 0.8 seconds when gets smaller than or equal to 10 -2 , which is when the D operator starts to be evaluated using a Taylor expansion. In this section we use our software Isolating_singularities to deal with a classical problem in robotics design and analysis. This problem is about determining the configurations in which a robot is allowed to move without damaging it; see for example However, regions that include the serial and parallel singularities (i.e., the boundaries of the generalized aspects) are left undecided since the numerical constrains are not satisfied. The results we present in section 8.3 can be considered as a continuation of the above-mentioned work in [CCG + 2012]. More precisely, we consider the manipulator RRRRR which we describe in Section 8.2. We certify the location of the geometric singularities (defined in Section 2.2.1) appearing in the reachable workspace and reachable joint-space (Definitions 8.1.1 (d)) as long as our strong assumptions are satisfied. Before that, we start with recalling preliminaries in Section 8.1.

Applications in robotic mechanisms

Preliminaries in robotics

Quoted from [LP2017, §2], a robot is mechanically constructed by connecting a set of bodies, called links, to each other using various types of joints. Actuators, deliver forces or torques that cause the robot's links to move. The positions that the robot attains are called configurations of the robot. The number of degrees of freedom of a robot is, roughly speaking, the smallest number of real-valued coordinates needed to represent the robot's configurations.

Robots can be classified into serial and parallel. Serial robots [Mer2006, §1.1] consist of series of links that are connected by joints equipped with motors. A parallel robot [Mer2006, §1.4.3] is made up of an end-effector with k degrees of freedom, and of a fixed base, linked together by at least two independent kinematic chains. Actuation takes place through k simple actuators.

Since the example we consider in our experiment is a parallel robot, we start with defining related vocabulary for modeling parallel robots: (c) [LP2017, Definition 2.1] For a robot with k degrees of freedom, the kdimensional variety Σ containing all possible configurations of the robot is called the configuration space (or C-space). The C-space Σ is embedded in the pose-command space. The C-space Σ of a robot is presented as a system of equations F (x, q) = 0 that characterizes the relation between the commands q (which represent the command variables) and the corresponding poses x.

(d) [CCG + 2012, §3.2] The projection of Σ on the pose space (resp. command space) is called the reachable workspace Σ x (resp. reachable joint-space Σ q ).

(e) [GA1990, Page 282] Let J F be the Jacobian matrix of F. Define the direct kinematic matrix A (resp. the inverse-kinematics B) to be the sub-matrix of J F obtained by removing the columns corresponding to q (resp. x).

(f) [GA1990, Page 282] The parallel (resp. serial) singularity set S parallel (resp. S serial ) is the set of points (x, q) in Σ such that det(A(x, q)) = 0 (resp.

det(B(x, q)) = 0). The serial (resp. parallel) singularities are also called first (resp. second) kinds of singularities.

(g) [CW1998, Definition 5] A generalized aspect A is defined as a maximal subset in Σ such that A is connected and for all (x 0 , q 0 ) ∈ A the matrices A(x 0 , q 0 ) and B(x 0 , q 0 ) are full-rank.

In this section, we assume that the dimensions of the pose and joint spaces are both equal to corank(J F (x, q)) for a generic (x, q) ∈ Σ.

Remark 8.1.1. It is important to emphasize that the notions of serial and parallel singularities differ from the singularity defined on Section 2.2.1. In the rest of this section, by geometric singularity we mean the one defined in Section 2.2.1.

To illustrate the objects that are defined above, we start with a simple example presented by Caro et al. [CCG + 2012]. Consider the planar robot PRRP depicted in Figure 8.1. A demo for simulating the motion of this robot is accessible in [Kra2021b] The robot consists of two prismatic joints sliding along two perpendicular axes. The positions of the prismatic joints are denoted by x (along x-axis) and q (along q-axis). These prismatic joints are connected through rigid links of lengths a, l and b respectively, with l is strictly larger than max(a, b). The three links are connected to each other by two revolute joints each of which is rotatable 360 degrees.

As illustrated in Figure 8.2, the pose-command space is hence R 2 and the configuration space Σ is the solution set of the equation (x -a) 2 + (q -b) 2 = l 2 . The pose space (resp. joint space) is [a -l, l + a] (resp. [b -l, l + b]). The serial singularity S serial is the set {(a + l, b), (a -l, b)}. Figure 8.4 illustrates the positions of the robot corresponding to the serial singularities. In these two configurations, the robot motion is restricted but not damaged. However, when the robot passes through a point in the parallel singularity S parallel = {(a, l+b), (a, b-l)} (the corresponding positions are depicted in Figure 8.3), a possibility of robot damage exists. Hence, in the design and the path planning processes of the robot, it is important to determine the location of the serial and parallel singularities in order to avoid any dangerous motion.

The green curve in Figure 8.2 is the union of all generalized aspects (Definition 8.1.1 (g)) of the robot. The robot is safely allowed to move from and to two points that are in the same generalized aspect (along that aspect). Notice that for every point (x 0 , q 0 ) in a generalized aspect, there exists a local bijection between the commands and the poses near that point. Equivalently, both of the matrices A(x 0 , q 0 ) and B(x 0 , q 0 ) are invertible. In the following section we consider a more challenging example, namely, RRRRR robot where the serial and parallel singularities are more complicated than in the above example.

Description of RRRRR robot

Consider the robot RRRRR (Figure 8.5) that consists of five links AP , P P, P P , P B and BA connected by five joints A, P , P, P and B. The links' lengths are respectively 8, 5, 8, 5 and 9. The angles BAP and ABP are the commands q 1 , q 2 respectively that vary in the interval (-π, π). The pose variables x 1 , x 2 are the coordinates of the point P and vary in the interval [-15, 15]. A simulation of this robot motion is accessible in [Kra2021a]. The kinematic equations are:    (x 1 -8 cos(q 1 )) 2 + (x 2 -8 sin(q 1 )) 2 -25 = 0, (x 1 -9 -5 cos(q 2 )) 2 + (x 2 -5 sin(q 2 )) 2 -64 = 0.

(8.1) Thus, the configuration space is a surface in the four-dimensional vector space represented by the variables x 1 , x 2 , q 1 and q 2 . Using the notation of Definitions 8.1.1 (e,f), the parallel (resp. serial) singularity is the set of points (x 1 , x 2 , q 1 , q 2 ) in the configuration space (i.e., satisfying (8.1)) that satisfy the equation det(A) = 0 (resp. det(B) = 0). In other words, the functions P 1 and P 2 that define the parallel and serial singularities are respectively:

P 1 =
[(x 1 -8 cos(q 1 )) 2 + (x 2 -8 sin(q 1 )) 2 -25, (x 1 -9 -5 cos(q 2 )) 2 + (x 2 -5 sin(q 2 )) 2 -64, det(A) = (2x 1 -16 cos(q 1 ))(2x 2 -10 sin(q 2 )) -(2x 2 -16 sin(q 1 ))

× (2x 1 -10 cos(q 2 ) -18)].

(8.2)

and

P 2 =
[(x 1 -8 cos(q 1 )) 2 + (x 2 -8 sin(q 1 )) 2 -25, (x 1 -9 -5 cos(q 2 )) 2 + (x 2 -5 sin(q2)) 2 -64, det(B) = (16(x 1 -8 cos(q 1 )) sin(q 1 ) -16(x 2 -8 sin(q 1 )) cos(q 1 ))

× (-10(x 2 -5 sin(q 2 )) cos(q 2 ) + 10(x 1 -5 cos(q 2 ) -9) sin(q 2 ))]

(

Results

As mentioned in the preamble of this chapter, our goal is to certify the nodes of the projections of the parallel and serial singularities with respect to the command and pose spaces. This is a continuation of classical work in robotics such as [CCG + 2012]. Applying the functions enclosing_curve and enclosing_singularities in software Isolating_singularities we get the following results:

The parallel singularity S parallel

All boxes that enclosing_curve returns are certified. Hence, P 1 satisfies Assumption A 1 . Taking the plane projections of S parallel with respect to (x 1 , x 2 ) and (q 1 , q 2 ) and solving Ball(P 1 ) accordingly we get:

• The plane projection of S parallel with respect to x 1 and x 2 :

The output of Isolating_singularities consists of:

(a) an empty set of uncertified solutions of Ball(P 1 ). Hence, P 1 satisfies Assumption ℵ 1 in B Ball , (a) The projection of S serial on R 3 with respect to (x 1 , x 2 , q 1 ). The projection of S serial on R 3 with respect to (x 1 , x 2 , q 2 ). (c) The projection of S serial on R 3 with respect to (q 1 , q 2 , x 1 ). with respect to (q 1 , q 2 , x 2 ). (f) The plane projection of S serial on R 3 with respect to (q 1 , q 2 ). (a) The projection of S parallel on R 3 with respect to (x 1 , x 2 , q 1 ). The projection of S parallel on R 3 with respect to (x 1 , x 2 , q 2 ). (c) The projection of S parallel on R 3 with respect to (q 1 , q 2 , x 1 ). (d) The projection of S parallel on R 3 with respect to (q 1 , q 2 , x 2 ). (f) The plane projection of S parallel on R 3 with respect to (q 1 , q 2 ). (b) an empty set of solution enclosures X = (x 1 , x 2 , y, r, t), with 0 ∈ t, and (c) a list of eleven Ball(P 1 ) solutions that are corresponding to nodes. As shown in Table 8.1, every row represents an enclosing box of a solution of Ball(P 1 ). Notice that the boxes returned by enclosing_singularities enclose only one solution of every pair of twin solutions. Moreover, recall that enclosing_singularities checks Assumptions ℵ 2 and ℵ 3 which is straightforward to verify by hand from Table 8.1. Hence, By Lemma 5.2.2, P 1 satisfies Assumptions A 1 , A 2 , A 3 , A 4 and A 5 .

Figure 8.8 shows the plane projection of S parallel with respect to x 1 and x 2 (brown boxes) in addition to the boxes enclosing the solutions of Ball(P 1 ) (black boxes). In fact, the black boxes are inflations of the output of Isolating_singularities since the width of the latter ones are not large enough to be visualized.

• The plane projection with respect to q 1 and q 2 :

In this case, enclosing_singularities returns an empty set of solutions (certified and uncertified) of Ball(P 1 ). Hence, Assumptions A 1 , A 2 , A 3 , A 4 and A 5 are satisfied. By Theorem 4.1.3, the plane projection of S parallel with respect to q 1 and q 2 is smooth; see Figure 8.9.

Table 8.1: The list of certified nodes enclosure that enclosing_singularities returns, where every row represents an enclosing box of a solution of Ball(P 1 ). Notice that enclosing_singularities returns exactly one solution of every pair of twin solutions.

q 1 q 2

Figure 8.9: The plane plotting of Isolating_singularities output for S parallel with respect to q 1 and q 2 .

The serial singularity S serial

Applying enclosing_curve to P 2 , both sets of certified and uncertified boxes are not empty (for eps min > 0). Hence, enclosing_curve does not certify the smoothness of S serial . Indeed, using a computer algebra system, we checked that the points

( 9 2 , √ 595
2 ) and ( 9 2 , -√

2 ) are both the plane projections (with respect to x 1 , x 2 ) of singular points of S serial . Projecting to the planes (x 1 , x 2 ) and (q 1 , q 2 ) and solving Ball(P 2 ) accordingly, we get:

• The plane projection with respect to x 1 and x 2 ; Figure 8.10: For all values of the parameter eps min that we chose, Isolating_singularities returns an empty set of certified Ball solutions. Moreover, the plane projections of uncertified boxes cover the projection of S serial . To understand this result, see Figure 8.11, where two different configurations have the same position P = (x 1 , x 2 ). This holds for uncountable number of configurations. Thus, Assumption A 2 is not satisfied (i.e., L n is uncountable). In fact, using a similar argument, we can see that L n is equal to S serial except the singular points.

• The plane projection with respect to q 1 and q 2 ; Figure 8.12:

In this case, enclosing_singularities returns an empty set uncertified Ball(P 2 ) solutions in addition to 14 boxes each of which encloses a certified node (Table 8.2). In Figure 8.12 we illustrate the projection of S serial with respect to q 1 and q 2 together with certified nodes which are inflated in the figure so they can be seen. The projection of unknown boxes (which contain geometric singularities of S serial ) are plotted in red. Table 8.2: The list of certified nodes enclosure that enclosing_singularities returns, where every row represents an enclosing box of a solution of Ball(P 2 ). Notice that enclosing_singularities returns exactly one solution of every pair of twin solutions.

q 1 q 1 q 1 q 2 q 2 q 2 Figure 8.12: An illustration of enclosing_singularities output for the serial singularity S serial of the robot RRRRR. The projection is with respect to the (q 1 , q 1 )plane (the command variables). The certification of the geometric singularities is done only where our strong assumptions are satisfied. The red boxes contain the projections of the (geometric) singular points of S serial .

Conclusion

We addressed in this thesis the problem of isolating the singularities of plane curves using certified and numerical methods. This is motivated by the fundamental problem of approximating (drawing) plane curves with the correct topology. We specifically considered the case of curves that are given as the plane projection of generic implicit curves in R n . The study of such curves was motivated by applications in robotics design. Our contribution in this thesis has two parts:

Theoretical contributions. Given a generic curve C contained in an n-dimensional closed box B, defined as the zero locus of a C ∞ -function from R n to R n-1 , we identified in Chapter 4 a regular square system, called the Ball system that encodes the singularities of the plane projection π C (C) of C. The idea of the Ball system is based on a theoretical analysis of generic curves in R n that classifies the points of C depending on the types of their projections. Thanks to its (square and regular) structure, this system is solvable using efficient certified numerical methods such as Interval Newton [Neu1991]. In Section 3.2, we precisely defined the conditions that C must satisfy in order to guarantee the correctness of our results. Furthermore, in the same section, we proved the genericity of these conditions via transversality theory. In the same way, we studied whether the silhouette of a generic surface in R n satisfies our conditions. In this case, we only proved the genericity of some of the assumptions and we stated the genericity of the missing part as Conjecture 3.2.18.

Algorithmic contributions and applications. We presented semi-algorithms that terminate if and only if the above-mentioned generic conditions are satisfied (Chapter 5). Furthermore, under the same assumptions, the semi-algorithms provide a set (or a superset) of isolating boxes of π C (C) singularities. We also presented an implementation of our semi-algorithms which we modified (by adding a threshold parameter eps min ) to ensure its termination even when our assumptions are not satisfied. Hence, given any curve C in R n , our software provides two sets of 2Dboxes: (i) a set of boxes each of which contains exactly one singularity, namely, a node of π C (C) and (ii) a set of boxes of widths smaller than eps min and above which our software is unable to verify the assumptions, and where, the behavior of π C (C) inside these boxes is thus unknown. Furthermore, any singularity of π C (C) is ensured to be in the union of these two sets of boxes even when our assumptions are not satisfied. We also discuss the performance of our software in Chapter 7 on non-algebraic instances for which no certified methods used to exist and algebraic instances that used to be out of reach by known certified methods. Finally, we finished this thesis with an application in robotics in Chapters 8.

Strengths, weaknesses and open problems

Our approach considers not only algebraic curves, but also those that are defined by C ∞ -functions as long as interval evaluations of the functions and their derivatives are provided. Moreover, our approach is the first certified numerical method that considers the curve C to be in R n for an arbitrarily large n. Our Ball-system-based approach is certified and appears to be efficient in several examples (see Chapter 7). Furthermore, it allows to isolate π C (C) singularities without computing an implicit equation of π C (C) in the plane for which the available approaches are usually restricted to the algebraic case and might increase the complexity of the method used.

Although, for technical reasons, we restricted C to live in a closed box in R n , we believe that our theoretical results in Chapters 3 & 4 still hold for any compact subset of R n . We also expect that the semi-algorithms in Chapter 5 can be developed to handle curves over compact sets as long as the membership problem17 over these sets is solvable. This generalization is, of course, helpful to tackle some visualization or robotics problems (that are similar to the ones that motivated us) where the region of interest is as general as compact sets.

It might also be interesting to investigate the possibility of generalizing the considered type of functions to the class C k , with k large enough, where C k is the class of functions that are differentiable k times with continuous derivatives. This generalizations is expected to produce the same types of singularities in the projection. In fact, our motivation behind this suggestion is purely mathematical.

Another interesting open problem is to investigate the possibility of computing a certified approximation of singular plane curves by first lifting the input curve to a smooth curve in R n , for some n, and then applying our approach. Even though it seems unintuitive to first lift for then to project, this approach is similar in spirit to blow-up techniques [JP2000, §5.3], although our lifting problem is global instead of local. However, if possible, it might allow to use certified numerical methods for computing the topology of implicit non-algebraic plane curves.

A similar idea of our approach has been already used by Diatta et al. [DMP2019]. The authors encoded the singularity of a surface in R 3 that is the projection of a generic surface in R 4 . This motivates the question of generalizing the idea of the Ball system to higher-dimensional objects. In particular, to find a regular system that encodes the singularities of the projection of a generic surface in R n to R 3 , with n > 4. If possible, this would help to address some problems raised when designing robots of degrees of freedom larger than two (see Section 1.1).

The genericity of the assumptions implies that a randomly chosen curve in R n satisfies our assumption. In this regards, it would be interesting to continue the investigation on whether the silhouette of a generic surface in R n satisfies our weak assumptions. We have already proved a part and we left the rest as a conjecture (Conjecture 3.2.18). Such silhouette curves are relevant because they naturally appear in parametric systems since they partition the parametric space with respect to the number of solutions of the system.

Recall that we added a parameter eps min in our semi-algorithm to enforce its termination. From a practical point of view, eps min also makes the algorithm practical with fixed-precision arithmetic since we can stop the refinements before we go beyond the limit of the fixed precision. In addition, it is useful if we want to ensure the termination even when not all assumptions are satisfied, certifying π C (C) singularities only on the part above which C verifies the assumptions. In particular, we can handle that way the projections of non-generic singular curves in R n .

The software we provided does not only offer singularity detection but also a visualization tool for π C (C). Figure 8.12 illustrates an example of how the output of our software look likes. Although the provided visualization does not guarantee the topology, the certified isolation of the singularities of π C (C) paves the way towards drawing π C (C) with the correct topology. Indeed, note that, for a given precision, if the drawing shows two branches that appear to intersect but without any reported box that contains the intersection, it is ensured that the two branches do not actually intersect; however, in order to determine the local topology (see Figure 9.1), further refinements are needed in our approach. We expect that an analogous approach to [IMP2018] might be found.

Our software handles only the functions that are implemented in the Minibex Boxes enclosing smooth parts of the projected curve A C B Figure 9.1: An example where our visualization tool does not provide a correct topology for some precision. A: Our software returns a set of boxes that encloses the projected curve without reporting any node (or unknown-status boxes). B & C: Two different scenarios of the plot in A after enough refinements. language18 (e.g., sin, cos, exp. . . ) which is the language used in our black-box solver. In fact, we also developed another version of our software where the user can implement interval evaluations of his/her own functions and their derivatives. However, this version is too slow in practice. An interesting project would be to investigate the possibilities of improving this version. For this goal, Python-FLINT [Joh2012b] offers an interesting toolbox for interval arithmetic.

Another natural open question regarding the semi-algorithms is its complexity. It is worth noticing that our semi-algorithms use a subdivision approach with diameter distance tests as defined by Burr et al. [BGT2020]. As such, it should be possible to study our complexities using the method of continuous amortization. This

1. 3

 3 Approximating the given curve preserving its topology, but loosing relative positions of the connected components. . . . . . . . . . . 1.4 An approximation of the plane curve that is ambient-isotopic to the given curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 An approximation of the plane curve where the vertices are boxes isolating singular or smooth critical points connected by segments approximating the curve (in blue). . . . . . . . . . . . . . . . . . 1.6 The curve defined by y 2 -(x 3 + x 2 ) = 0. . . . . . . . . . . . . . 1.7 Illustration of a node and cusps in the plane projection of a smooth curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Left: An example of a box in R n that satisfies (a) and whose projection contains a node. Right: A simplified illustration of pruning the domain of the Ball system. The brown and orange boxes stratify (a), also the pair of red boxes satisfy (b). Using these boxes, the blue boxes in R 2n-1 are computed. . . . . . . . . . . . . . . . . . 2.1 An example in R 2 where f and g define the same germ. . . . . . . xiii 2.2 Left: At an A 1 singularity, two branches of the curve intersect transversally. Right: At an A 2k+1 singularity with k > 1, the tangent lines of the two branches at the intersection point coincide. . 3.1 Illustration of the assumptions. . . . . . . . . . . . . . . . . . . . 3.2 Illustration of the proof of Lemma 3.2.6. . . . . . . . . . . . . . . 3.3 Proof of Lemma 3.2.7: The first case where v and v are independent.

4. 1

 1 Illustration of a node and cusps in the plane projection of a smooth curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The curve C (red) and its plane projection π C (C) (blue) of Example 4.1.11 displaying a cusp singularity. . . . . . . . . . . . . . . 4.3 The curve C (red) and its plane projection π C (C) (blue) of Example 4.1.12 displaying a node singularity. . . . . . . . . . . . . . . 4.4 An example of Lemma 4.3.1 for n = 3. . . . . . . . . . . . . . . 4.5 Comparing the singularities considered in Lemma 4.3.4 with A 2k+1 types. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5. 3

 3 An example of B and B satisfying (2) in Lemma 5.3.1. The projection has a loop. This implies that there are a x 1 -critical point and a x 2 -critical point in the union of B and B . In this example, these points are different and the projection has a node. . . . . . .5.4 A simplified illustration of pruningB Ball as described in Proposition 5.3.3. The red boxes satisfy (1). On the other hand, the brown (resp. orange) box satisfies (2). Hence, we solve the Ball system over f Ball pre-images (in blue) of the red boxes, the brown box and the orange box. . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Experiment 1: Plane projection of an analytic curve in R 3 with one node and one ordinary cusp. . . . . . . . . . . . . . . . . . . . . 7.2 Experiment 2: Plane projection of an analytical curve C in R 4 . Each of the 43 black boxes contains a node of π C (C) and is the projection of a box in R 7 containing one solution of Ball(P ). . . . 7.3 Experiment 3: High degree algebraic curve in R 3 generating 7 nodes.115 8.1 A non-singular configuration of PRRP where there is a local bijection between poses and commands. . . . . . . . . . . . . . . . . . 8.2 Geometric modelling of PRRP motion. . . . . . . . . . . . . . . 8.3 The parallel singular poses of PRRP where the pose x admits two commands. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4 The serial singular poses of PRRP where the command q admits two poses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5 RRRRR robot in a non-singular configuration. . . . . . . . . . . . 8.6 Illustration of the projections of S serial on different variable subspaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7 Illustration of the projections of S parallel on different variable subspaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.8 Enclosing the plane projection of S parallel with respect to x 1 and x 2 .

Figure 1

 1 Figure 1.1: Left: Illustration of a configuration space of a 2-degrees-of-freedom robot. In this example, it is a surface embedded in R 3 . Generically, the silhouette of this surface is a (smooth) curve. Right: Projection of the silhouette on the motor space.

  Figure1.2: An approximation of the curve in Figure1.1 (right) that preserves the topology and the locations of the critical points. The approach used for this goal is to isolate the critical points, and then, to connect them properly using the stateof-the-art methods since the remaining part is smooth. Additional black vertices are added to improve the visualization. The number on each face represents the number of solutions of the direct kinematic problem.

  Figure 1.5: An approximation of the plane curve where the vertices are boxes isolating singular or smooth critical points connected by segments approximating the curve (in blue).

Figure 1

 1 Figure 1.6: The curve defined by y 2 -(x 3 + x 2 ) = 0.

Figure 1 . 7 :

 17 Figure 1.7: Illustration of a node and cusps in the plane projection of a smooth curve.

Figure 1

 1 Figure1.8: Left: An example of a box in R n that satisfies (a) and whose projection contains a node. Right: A simplified illustration of pruning the domain of the Ball system. The brown and orange boxes stratify (a), also the pair of red boxes satisfy (b). Using these boxes, the blue boxes in R 2n-1 are computed.

Figure 2

 2 Figure 2.2: Left: At an A 1 singularity, two branches of the curve intersect transversally. Right: At an A 2k+1 singularity with k > 1, the tangent lines of the two branches at the intersection point coincide.

  Having thus defined convergent sequences, we now know what are the closed and open sets in

  and assumptions . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1.2 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 Genericity of the assumptions . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.1 Genericity of the assumptions for a curve in R n . . . . . . . . . . 41 3.2.2 Genericity of the assumptions for the silhouette of a surface in R n . 49

  For a positive integer n, a closed (resp. an open) n-box is the Cartesian product of n closed (resp. open) interval. Assume that n 3 and let B be an open n-box and B be the topological closure of B with respect to the usual topology

  (b) Assumptions A2 and A4: the sets Lc and Ln are finite and do not intersect the boundary of B. Assumption A3: No three points (counted with multiplicity) have the same projection. Assumption A3: No three points (counted with multiplicity) have the same projection. Assumption A5: points in πC(C) are only smooth or nodes.(e) Assumption A -5 : points in πC(C) are either smooth, nodes or ordinary cusps.

Figure 3 . 1 :

 31 Figure 3.1: Illustration of the assumptions.

Figure 3

 3 Figure 3.2: Illustration of the proof of Lemma 3.2.6.

Figure 3

 3 Figure 3.3: Proof of Lemma 3.2.7: The first case where v and v are independent.

Figure 3

 3 Figure 3.4: Proof of Lemma 4.3.3: The second case where v and v are co-linear.

Figure 3

 3 Figure 3.5: Example: The silhouette of the torus in R 3 with the coordinates (x, y, z) (source: [IMP2016a]) .

Contents 4. 1

 1 Encoding the singular points of the plane projection . . . . . . . . . . . 57 4.2 Singularities induced by Ln . . . . . . . . . . . . . . . . . . . . . . . . 64 4.3 Singularities induced by Lc . . . . . . . . . . . . . . . . . . . . . . . . 64 4.4 Regularity of the Ball system . . . . . . . . . . . . . . . . . . . . . . . 73

Figure 4 . 1 :

 41 Figure 4.1: Illustration of a node and cusps in the plane projection of a smooth curve.

Figure 4 . 2 :

 42 Figure 4.2: The curve C (red) and its plane projection π C (C) (blue) of Example 4.1.11 displaying a cusp singularity.

Figure 4 . 3 :

 43 Figure 4.3: The curve C (red) and its plane projection π C (C) (blue) of Example 4.1.12 displaying a node singularity.

Figure 4 . 4 :

 44 Figure 4.4: An example of Lemma 4.3.1 for n = 3.

Figure 4 . 5 :

 45 Figure 4.5: Comparing the singularities considered in Lemma 4.3.4 with A 2k+1 types.

Contents 5 . 1

 51 Interval arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.2 Semi-algorithms to check assumptions and isolate singularities . . . . 85 5.3 Semi-algorithms improvements . . . . . . . . . . . . . . . . . . . . . . 92 5.3.1 Subdivision-based solver on a high-dimensional system . . . . . . 92 5.3.2 Evaluating the operator D . . . . . . . . . . . . . . . . . . . . . . 98 5.3.3 Improvement integration . . . . . . . . . . . . . . . . . . . . . . 99

Semi-algorithm 3 :

 3 Checking the weak assumptions and computing a superset of the singularities of π C (C) Input: An open n-box B and a smooth function P from B to R n-1 . Termination: If and only if P satisfies Assumptions A 1 , A 2 , A 3 , A 4 and A - 5 in B.

  A 1-5 and isolating the singularities of π C (C) Input: An open n-box B and a smooth function P from B to R n-1 .

  (a) every box B in enclosing_curve contains points of C, (b) the intersection of two distinct B, B in enclosing_curve lives in their boundaries, moreover, the set B \ B still contains points of C, (c) for non-disjoint B and B in enclosing_curve, the curve C is connected in B ∪ B , then every node in π C (C) is contained in at least one of the following 2-boxes:

Figures 5 .

 5 Figures 5.2 and 5.3 show two examples of (2) where B = B .

Figure 5 . 2 :

 52 Figure 5.2: An example of B and B satisfying (2) in Lemma 5.3.1. Their projections have a node in the boundary. The curve C is tangent to the boundaries of bothB and B at q and q respectively. The curve and the boxes intersect at q such that q does not have the same plane projection as q and q .

Figure 5 . 3 :

 53 Figure 5.3: An example of B and B satisfying (2) in Lemma 5.3.1. The projection has a loop. This implies that there are a x 1 -critical point and a x 2 -critical point in the union of B and B . In this example, these points are different and the projection has a node.

  2 and t = [0, max( By-By 2 4 )]}. (b) There exist (possibly equal) boxes B = B (x 1 ,x 2 ) ×B y and B = B (x 1 ,x 2 ) ×B y in enclosing_curve satisfying (1) or

Figure 5

 5 Figure5.4: A simplified illustration of pruning B Ball as described in Proposition 5.3.3. The red boxes satisfy (1). On the other hand, the brown (resp. orange) box satisfies (2). Hence, we solve the Ball system over f Ball pre-images (in blue) of the red boxes, the brown box and the orange box.

1 2 ( 4 .

 24 B y + B y ), r = [-1, 1] n-2 if B ∩ B = ∅, otherwise, r = By-B yBy-B y , t = By-B y 2 The arithmetic operations mentioned in Proposition 5.3.3 are the natural box function of the usual arithmetic operations defined on R n (see Example 5.1.2). Figure 5.4 illustrates a global picture of the steps above.

1:

  Compute a set enclosing_curve of closed n-boxes that enclose C (using a black-box solver) with Properties (a), (b) and (c). 2: Find the boxes of enclosing_curve satisfying (1) or (2) (Remark 5.3.2) 3: Solving_domains := {} 4: for a box B = B in enclosing_curve satisfying (2) do 5: (x 1 , x 2 , y) = B, r = [-1, 1] n-2 and t = [0, max( By-By 2 4

  -party libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.2 The Python software Isolating_singularities . . . . . . . . . . . . . . . 107

  4673, 2.4675] U = {Box2 = [-0.06, 0.04] × [-0.04, 0.07] × [-3.15, -3.14] × [1, 1] × [-0.01, 0.01]}.

Figure 7 .

 7 Figure 7.2 illustrates the output of our improved Semi-algorithm 4 for the curve defined by

  Figure 7.2: Experiment 2: Plane projection of an analytical curve C in R 4 . Each of the 43 black boxes contains a node of π C (C) and is the projection of a box in R 7 containing one solution of Ball(P ).

Figure 7 .

 7 Figure 7.3 illustrates the 7 nodes of the projection of C and running times are given in Table 7.1.Note that since P is polynomial, the implicit equation of π C (C) can be computed

  Table 7.2: Experiment 4: Performances for different values of .

  Contents 8.1 Preliminaries in robotics . . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.2 Description of RRRRR robot . . . . . . . . . . . . . . . . . . . . . . . . 122 8.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 8.3.1 The parallel singularity S parallel . . . . . . . . . . . . . . . . . . 123 8.3.2 The serial singularity S serial . . . . . . . . . . . . . . . . . . . . 128

  [BL1986,Cha2010,CCG + 2012]. Caro et al. [CCG + 2012] provided an approach, based on branch and prune algorithm, to compute a certified enclosure of the generalized aspects (see Definition 8.1.1 (g)) which fully automatically applies to arbitrary robot kinematic model under numerical constrains [CCG + 2012, §4.1].

  Definitions 8.1.1. Consider the following notions. (a) ([Mer2006, §1.4.3]) A fully parallel robot is a parallel robot for which the number of chains is equal to the number of degrees of freedom of the endeffector. (b) [CCG + 2012, §3.2] The pose-command space is the product of the pose and command spaces which are ambient spaces of all poses and all commands respectively.

Figure 8

 8 Figure 8.1: A non-singular configuration of PRRP where there is a local bijection between poses and commands.

Figure 8

 8 Figure 8.2: Geometric modelling of PRRP motion.

Figure 8

 8 Figure 8.3: The parallel singular poses of PRRP where the pose x admits two commands.

Figure 8 . 4 :

 84 Figure 8.4: The serial singular poses of PRRP where the command q admits two poses.

Figure 8 . 5 :

 85 Figure 8.5: RRRRR robot in a non-singular configuration.

  The projection of S serial on R 3

  The plane projection of S serial on R 3 with respect to (x 1 , x 2 ).

Figure 8 . 6 :

 86 Figure 8.6: Illustration of the projections of S serial on different variable sub-spaces.

  The plane projection of S parallel on R 3 with respect to (x 1 , x 2 ).

Figure 8 . 7 :

 87 Figure 8.7: Illustration of the projections of S parallel on different variable subspaces.

Figure 8

 8 Figure 8.8: Enclosing the plane projection of S parallel with respect to x 1 and x 2 .

Figure 8

 8 Figure 8.10: Plotting the result of applying Isolating_singularities to S serial with respect to x 1 , x 2 , where it plots only red boxes since neither Assumption A 2 nor Assumption A 1 is satisfied.

Figure 8

 8 Figure 8.11: An example of two different (serial) singular configurations of RRRRR that define the same pose (x 1 , x 2 ).

  Experiment 4: Performances for different values of . . . . . . . .

	Contents
	List of Tables
	7.1 Timings (in seconds) and numbers of boxes in Experiments 1 to 3.
	7.2
	Bibliography
	x	ix

8.1 The list of certified nodes enclosure that enclosing_singularities returns, where every row represents an enclosing box of a solution of Ball(P 1 ). Notice that enclosing_singularities returns exactly one solution of every pair of twin solutions. . . . . . . . . . . . . . . 8.2 The list of certified nodes enclosure that enclosing_singularities returns, where every row represents an enclosing box of a solution of Ball(P 2 ). Notice that enclosing_singularities returns exactly one solution of every pair of twin solutions. . . . . . . . . . . . . . . xi

  where Reg(W ) (resp. Sing(W)) is the set of smooth (resp. singular) points in W. Recall that codim(Sing(W )) > n since Sing(W ) is a proper closed sub-variety of W [BCR1998, Proposition 3.3.14]. By Corollary 2.2.16, there exists a residual set Γ

  If Assumptions A 1 , A 2 , A 3 and A 4 are satisfied, then Assumption A - 5 is satisfied if and only if Assumption A - 5 is satisfied for all

			This implies
	that A is equivalent to	0 1 1 0	and hence is invertible.
	Lemma 4.4.7.

  Assumption A 3 is satisfied. Now, Since Assumptions A 1 ,A 2 , A 3 and A 4 are satisfied and since all solutions of Ball(P ) are regular, by Theorem 4.4.2, we have that Assumption A - , we are ready to check Assumptions A 2 , A 3 , A 4 and A - 5 using ℵ 1 , ℵ 2 and ℵ 3 . Since Lemma 5.2.2 requires Assumption A 1 , we start by checking that assumption with Semi-algorithm 1 that is based on subdivision.

	5
	is also satisfied.
	Using Lemma 5.2.2Semi-algorithm 1: Checking Assumption A 1

Input: An open n-box B and a function P from B to R n-1 . Termination: If and only if P satisfies Assumption A 1 in B. Output: True 1: L := {B} 2: while L = ∅ do 3:

  . Recall that every box U of L (after the change in Item i) is computed in such a way that f Ball (U) is a pair of boxes in enclosing_curve. Hence, no need to check the condition ( f Ball (U))∩(B×B) = ∅ in Step (5) of Semi-algorithm 2 since it is never satisfied. iv. As mentioned in Section 5.3.2, when solving the Ball system (Semialgorithm 2, Step (8)) over a box -inflation(U), we consider the two different approaches to evaluate the operator D. More precisely, if the tinterval [t min , t max ] of -inflation(U) contains zero, 10 we use Taylor form over the sub-box of -inflation(U) satisfying t min t min(δ Taylor , t max ). If t max > δ Taylor , we solve the Ball system over the complimentary sub-box (i.e., with δ Taylor t t max ) as defined in (4.1).

  6.1.1. Detecting pairs of boxes in enclosing_curve satisfying (1) or (2) in Lemma 5.3.1 can be done using a black-box tool. It is a software that reports the pairs of intersected boxes of a given set of boxes. In our implementation, we use the SciPy class KDTree [VGO+ 2020].15 This class provides a tool for quick nearest-neighbor lookup for points in R n . Considering the center of every box in enclosing_curve, we can find the centers of the other boxes that intersect it. Although KDTrees tools offer a fast way to detect the neighbored boxes, it might return in addition spurious boxes which we simply remove by checking whether a pair of boxes returned by KDTree intersect.

In addition to the previously-mentioned tools, other libraries are used such as SymPy [MSP + 2017], NumPy [HMvdW + 2020] and Matplotlib

[Hea2007]

.

  Table 7.1: Timings (in seconds) and numbers of boxes in Experiments 1 to 3.

			Enclosing curve C	Solving Ball	Output
						system
		Boxes
	Experiment 1 Experiment 2 Experiment 3 Experiments	0.03 0.01 0.1 0.03 0.01 0.1 0.03 0.01 max. width 0.1	1835 5427 2243 6759 19583 1151 2503 6347 Tree size 535	Output boxes 456 1354 520 1639 4847 203 523 1482 134	Time 0.3 0.7 1.1 3.4 10.2 1.0 1.8 4.5 0.1	Tree size 90 188 6098 1078 372 655 272 163 70	3.8 4.4 52 35.4 4.2 3.5 5.7 35.8 Time 3.6	5.3 Total time 5.2 10.2 4.1 5.1 53.1 38.8 45.6 3.7	Singularity 2 43 7 boxes

  8.3) Figures 8.6 and 8.7 illustrate the projections of the parallel and serial singularities on different variable sub-spaces. These figures are plots of the boxes returned by enclosing_curve.

Topology. We first consider "simple" properties such as the number of connected components, intersection points, etc. More formally, we wish to find a plane graph that is homeomorphic to the given curve. Figure1.3 (left) shows a plane curve to be approximated. An approximation of the latter curve that fulfills these conditions

Wolfgang Gröbner was Buchberger's thesis advisor.

http://www.ibex-lib.org/doc/minibex.html.

Recall that 2n -1 is the size of the Ball system.

The software is accessible at https://github.com/gkrait/Isolating_singularities.

Note that the converse is not true as the vertical (double) line defined by x 2 1 = x 2 = 0 in R 3 is smooth but the rank of its Jacobian is never full.

√ t (P i (q 1 ) -P i (q 2 )) = 0, for all i ∈ {1, . . . , n -1}. Thus, X is a solution of Ball(P ).

For a box B and > 0, -inflation(B) is the box that has the same center as B and its width is that of B multiplied by (1 + ). The box -inflation(B) thus contains B, the exact value of is not important for the algorithm and it is usually set to 0.1 in subdivision algorithms[Rum2010].

Regarding the following properties, Items (a) and (b) are to remove redundancy which, of course, improves the perfor-

Notice that the x 1 -values of q, q and q might be equal. This case occurs when C is a segment orthogonal to R 2 in B ∪ B .

In fact, a black-box tool in our implementation uses fixed precision arithmetic; see Section 6.1 for more details.

Notice that, in this case, t min is in (-∞, 0] by -inflation.

The software is accessible at https://github.com/gkrait/Isolating_singularities.

The software is accessible at https://github.com/gkrait/Isolating_singularities.

That is, when our software succeeds to verify the weak assumptions (Section 3.1.2).

Recall that eps min is the parameter used in Semi-algorithm 6.

This reference is related to SciPy which is a collection of open source software for scientific computing in Python where KDTree is a class of this collection; see https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html.

The software is accessible at https://github.com/gkrait/Isolating_singularities.

The membership problem asks whether a given point is in a given subset in R n .

http://www.ibex-lib.org/doc/minibex.html.

Remerciements

Part II

Software, experiments and application in robotics Conclusion

should yield explicit bounds for the case of polynomial input either in the worst case (as in [BGT2020]), or in a smoothed complexity setting (as in [CETC2019]).