{'fraction_non_alphanumeric': 0.04561166514362907, 'fraction_numerical': 0.010157998596156648, 'mean_word_length': 3.760307653511537, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 46, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Robust and Efficient Method for Solving Geometrical Constraint Problems by Homotopy

The goal of Geometric Constraint Solving is to find 2D or 3D placements of some geometric primitives fulfilling some geometric constraints. The common guideline is to solve them by a numerical iterative method (e.g. Newton-Raphson method). A sole solution is obtained whereas many exist. But the number of solutions can be exponential and methods should provide solutions close to a sketch drawn by the user. Assuming that a decomposition-recombination planner is used, we consider irreducible problems. Geometric reasoning can help to simplify the underlying system of equations by changing a few equations and triangularizing it. This triangularization is a geometric construction of solutions, called construction plan. We aim at finding several solutions close to the sketch on a one-dimensional path defined by a global parameter-homotopy using the construction plan. Some numerical instabilities may be encountered due to specific geometric configurations. We address this problem by changing on-the-fly the construction plan. Numerical results show that this hybrid method is efficient and robust.

Geometric Constraints Solving Problems (GCSP) arise in many fields such as CAD, robotics or molecular modeling. This is to determine the positions of geometric elements (points, lines, planes, circles, etc.) that must satisfy a set of constraints such as distances, angles, tangencies and so on. Commercial solvers generally rely on numerical methods (Newton or quasi-Newton) that provide a single solution. But even when problems are well-constrained, there may be several interesting solutions. Example of Fig. 1 shows a 2D sketch drawn by a user who imposed five constraints on it: 4 distances and 1 angle. For some assignments, there could be two solutions very close to each other (Fig. 1b) and there is no reason for excluding one from another. However, providing all the solutions is not always relevant since the number of solutions can be significant. In a 2D triangle "mesh" with n well-constrained triangles, a complete solving leads to 2 n solutions, some of which are very far from the shape of the sketch. Our aim is to design a method that uses the sketch to guide the research of several solutions. We assume that the considered GCSP are structurally well-constrained. Roughly speaking, this means that there exists some assignments for dimensions leading to finitely many solutions. Several methods can yield more than one and even all the solutions. Subdivision methods [15,5] give good results but the number of boxes to be explored can be huge. In algebraic approaches, homotopy methods have been successfully studied in this area [3] but only for small size problems. Indeed, the number of homotopy paths to follow grows exponentially with the number of constraints. [11] proposes to use the sketch to define a parameter-homotopy. Thus a sole path is followed but a sole solution is obtained.
Another way to get several solutions comes from geometric methods. Using geometric construction rules a construction plan is produced by the solver. The construction plan is evaluated numerically to yield different solutions. However, some problems in 2D and most 3D problems can not be solved by a sequence of elementary geometric constructions. To circumvent this, [6] proposes a reparameterization: constraints are replaced by some constraints of distances or angles in order to make possible the production of a symbolic construction plan. Next, the numerical resolution must find the values of added constraints to satisfy the deleted constraints. Again, the choice of the numerical resolution method is crucial to provide several solutions. [8] presents a first attempt for using a homotopy method along with a parameterized construction plan. The idea is to follow a homotopy path to which belongs the sketch.
In [9] it is shown that this path may contain several solutions and that parameterized construction plans can be used to detect solutions at infinity. In addition, the construction plan can be used to generate new sketches and thus new paths.
By combining homotopy and reparameterization, the number of unknowns of the problem is greatly reduced and these unknowns apply only to added constraints. The solution is structured in a construction plan, and the price to pay is that this plan is much longer to evaluate than polynomials of constraints. This method provides very good results, both in execution time and number of solutions found, but is unstable at the boundary of the domain. The boundary corresponds to geometric situations which are clearly identified such as a circle-circle intersection where the circles are almost tangent. In these situations, numerical path tracking is difficult. This paper justifies the resolution by homotopy of a reparameterized construction plan. It defines the boundary of the domain of definition and characterizes it at specific geometric configurations. The originality of this work is a change of the construction plan on-the-fly to avoid getting too close to the boundary. In this way we solve the stability problems of this method.
This paper is structured as follows. Section 2 reviews homotopy continuation methods and outlines our method on an example. Section 3 introduces definitions and basic concepts on construction plans. Section 4 presents homotopy paths tracking on construction plan. Section 5 explains how to change a construction plan on-the-fly to avoid critical situations. Section 6 presents some results and section 7 concludes this work.
2 Homotopy solving 2.1 A brief on classical homotopy methods Homotopy methods use continuation approach to find zero values of F : K m → K m by embedding it in an homotopy function H that interpolates an initial function F 0 : K m → K m into F . H can for instance be defined as H(x, t) = (1 − t)F (x) + tF 0 (x). When H is assumed to be regular and defined on K m × [0, 1], connected components of the set H −1 (0) = {(x, t) ∈ K m ×[0, 1]|H(x, t) = 0} are one dimensional manifolds which borders are elements of hyperplanes P i = {(x, t) ∈ K m × R|t = i} for i ∈ {0, 1}. These connected components are called homotopy paths and are diffeomorphic to a close, a semi-open or an open interval, or a circle.
If a zero value x 0 of F 0 is known, a homotopy paths S is followed from (x 0 , 0) by a path tracking method, until a zero-value of F is found provided that (a) S is diffeomorphic to [0,1] and has its other extremity in P 1 . Bad cases that can occur in a general continuation framework are the following: (b) S is diffeomorphic to [0, 1] but has its other extremity in P 0 and no zero values of F is found, or (c) S is diffeomorphic to [0, 1[, has infinite length and the tracking process does not terminate ( Fig. 2 shows these different cases).
Most of literature concerning homotopy methods addresses situations where K = C, and components of F are polynomials; see [17] for a survey. We introduce here some selected points.
When K = C is assumed, it is possible to ensure with probability one that H is regular by using the so-called gamma-trick and that homotopy paths of H are strictly increasing with respect to the interpolation parameter t. Hence case (b) can not occur.
When components of F are polynomials one can construct F 0 having at least as many roots as Inria Unknowns: point P0, ..., P5 Parameters: length h0, ..., h8 Constraints: distance(P0, P1) = h0 distance(P1, P2) = h1 ... distance(P1, P4) = h7 distance(P2, P5) = h8 00 00 11 11 0 0 1 1  F . In this case all roots of F are linked to roots of F 0 , and paths (c) cross P 1 in a point of infinite norm, called solution at infinity. Such paths can be detected before following them in a complex projective space by using an homogenization of H. Finding all complex roots of a polynomial map F by homotopy is then conditioned to the ability of constructing an initial polynomial map F 0 admitting as many known roots as F . The number of roots of F is bounded by the total degree of F , but often over-estimated from far by the latter, especially when polynomials are sparse and structured. Starting an homotopy with this initial system leads to follow too much path. The bound of the number of roots can be sharpened by exploiting an homogeneous structure of F (see [12]) or thanks to the Bernstein's theorem [7]. A survey on path tracking processes can be found in [1]. The most common one uses a prediction correction approach: knowing an initial (approximated) point (x 0 , t 0 ) on S, a predictor compute a new approximated point of S, and a corrector brings it much closer to the path. Prediction is often done along the tangent T on (x 0 , t 0 ) of S with a prediction step δ, and the corrector consists of Newton-Raphson iterations on H that has to be squared thanks to an additional equation. The latter can constraint the correction to lie on the hyperplane orthogonal to T . When F is polynomial and S is proved to be monotonic with respect to t, corrected point is searched such that t = t 0 +δ. The step δ is classically adapted dynamically trough the tracking process. Initialized by [10], recent researches on path tracking attempt to improve robustness and goes trough certification by using clever arithmetic such multi-precision arithmetic ( [2]) or interval arithmetic ( [10,4,13]).

Translating a GCSP into a system of equations (F)
where C is a set of constraints, X a set of unknowns and A a set of parameters. Usually, unknowns are geometric entities such as points and lines and parameters are dimensions such as lengths or angles. Constraints are predicative terms such that distance(p1, p2) = k for a distance constraint between points p1 and p2. A example is given in Fig. 3 where the problem consists in constructing six points P 0 , . . . , P 5 in a plane knowing 9 distances pairwise. We call this GCSP K 3, 3 1 . This definition of GCSP is general and for some statements it could be useful to consider some geometric entities as parameters and some dimensions as unknowns.
Considering that a numerical function is associated to each constraint, a GSCP is easily translated into a numerical function where algebraic unknowns are associated to GCSP unknowns. Since point-point distances constraints are invariant up to rigid motions, solutions are sought in a chosen reference; for instance in K 3,3 we could search figures with P 0 at the origin and P 1 with null ordinate. Algebraic unknowns {x 0 , . . . x 8 } are assigned to remaining free coordinates of GCSP unknowns: x 0 to the abscissa of P 1 , (x 1 , x 2 ) to the coordinates of P 2 , and so on.
Then, equation F 0 (x 1 ) = P 0 P 1 − h 0 where P 0 P 1 denotes euclidean distance between P 0 and P 1 is associated to distance constraint distance(P 0 , P 1 ) = h 0 . We will note F : R 9 → R 9 the function whose components are numerical functions associated to constraints of G, and call F the numerical function associated to G. Since G is structurally well constrained, F is well-formed: it has as many variables as equations, and no subset of equations contains fewer unknowns than equations.
In a numerical function F , we also denote by X the set {x 0 , . . . x m } of algebraic unknowns, and zeros of F are called figures. Since F also depends on parameters A, and we make it explicitly appear by noting If A so are the values of parameters given by the user, solving G for parameters A so is equivalent to solve the system of unknowns X: In homotopy path tracking, a system whose some solutions are known is required. As previously mentioned we will use the sketch as initial solution. Values for X read on the sketch drawn will be referred to as X sk . We have especially where A sk is the values of parameters read on the sketch.

Searching zero-values of F by homotopy
When components of F are polynomials, all solutions in C m of equation (2) can be provided by a classical homotopy method. [3] proposes to first use symbolic reduction on F to obtain F : C m → C m with m < m in order to decrease the number of path to be followed. But most of found solutions are complex whereas in applications only real figures are relevant. GCSP generally admit a number of solutions that grows exponentially with the number of constraints; in CAD applications, a user may prefer obtaining one or several solutions close to the sketch he drawn (that is supposed to represent well the desired solution) than choosing between an exponential number of figure. With this motivation, [11] uses the sketch to build the homotopy function H: taking (3) into account, H is defined as H(X, t) = F (X, tA so +(1−t)A sk ), and (X sk , 0) is the sole known solution to H(X, t) = 0. Only one homotopy path is followed, in the real space R m × [0, 1], and only one solution is found, but it is considered to be close to the sketch, as a continuous deformation of the latter.
The main limitations of a real homotopy lies in the existence of homotopy paths (b) and (c): following only one path may lead to obtain a solution of H(X, 0) = 0 in a finite time, or no solution in an infinite time. In the previous paper [9], we have proposed to re-use the main idea of [11], that is building the homotopy function from the sketch and tracking paths in a real space in order to obtain only real solutions close to the sketch. We have showed how to define the homotopy function in such a way that paths (c) can be detected during the tracking process because figures lying on it converge to a special geometric configuration. Moreover others paths, Inria considered in R m × R, are showed to be diffeomorphic to circles, and can be tracked until they loop, to obtain more solutions. Finally it shows how to use constructive geometry to find more solutions.
If this approach does not ensure to obtain neither one nor all the solutions of a GCSP, it is efficient in practice and sometimes allows to get all solutions. It provides several solutions to problem that are too big to be solved with a classical homotopy method. For the sake of selfsufficiency we recall here main contributions of [9] on which relies the method presented here. The idea of tracking a path in R m × R to obtain several solutions on a sole path is not new, and is called global homotopy in page 81 of [1].

Using a global homotopy to solve GCSP
Key results of our global homotopy approach comes from the way the parameters A sk are interpolated into A so . Let us introduce two notions in order to characterize it.
With GCSP K 3,3 of figure 3, with A sk the sketch values and A so the required values, the function a : R → R 9 is made up of components a i that could be for instance: with c > 0 is an interpolation function from A sk to A so .

Definition 2
If a is an interpolation function between α ∈ R and β ∈ R, we call positive support of a, and note it Supp(a), the set {x ∈ R|a(x) ≥ 0}.
The map a 0 defined in equation (4) has a positive support that is equal to the close real interval between the two real roots of a 0 .
If a is an interpolation function from A sk to A so and F a numerical function associated to a GCSP, we define the homotopy function H : R m × R → R m as: H(X, t) = F (X, a(t)).
For instance, the homotopy function associated to the GCSP K 3,3 is H : where a is defined as in equation (4). Main contributions of [9] are resumed by the following affirmation. If the following assumptions are met: (i) H is well defined over a m + 1 dimensional manifold DH ⊆ R m × R without any border,  Figure 4 presents such different paths. If an homotopy path S is diffeomorphic to an open interval then extremities of S converge either to a point in (R m × R) \ DH either to a solution at infinity. These cases correspond to special geometric configuration that can be detected during the tracking process. Assumption (i) holds when constraints are distance, colinearity, ..., angles between three points. The numerical function associated to the latter is not defined when two of the three points are coincident, and this geometric configuration can characterize a bounded limit of a homotopy path diffeomorphic to an open interval. Assumption (ii) stands for almost all values of A sk and A so , when each constraint depends on a numerical parameter. We also suppose that parameters are independents, that is each parameter is involved in a sole constraint. Assumption (iii) means that G involves at least one distance constraint, this is stronger than necessary but ensures the result to be true.
This gives rise to the simple following procedure, which is parameterized by a path tracking method: Continuation solving of GCSP Input: GCSP C[X, A], values A so , sketch X sk Output: list of solutions L sol 1. Compute values A sk from X sk , build an interpolation function a from A sk to A so , build H.
2. Follow path of H −1 (0) from (X sk , 0) until it loops on (X sk , 0) while checking if : • it converges to a solution at infinity, then stop; • it converges to a figure outside of DH, then stop; • it crosses hyperplane P 1 , then append solution to L sol .
The amount of computation required by the step 2 to follow a path S depends of course on the chosen prediction step and on the arc-length, but is also linked to the cost of an iteration. If the path tracker that is used is a simple predictor-corrector with a first-order predictor and a correction performed by Newton-Raphson iterations 2 , it requires at each iteration several inversions of an m × m matrix. The cost of our global homotopy solving is then directly related to the number of components of H.
The point of the present paper is to use constructive geometry to track paths of H in a space of dimension d, with d m.
2.4 Using reparameterization to quasi-triangularize the system of equations The presentation is made with example K 3,3 , but what is said is directly generalizable to others GCSP.
Consider the function H R defined as follows: where H is defined as in equation (6) and let us state the two following remarks.
Remark 1 Let k sk be the length P 0 P 2 read on the sketch X sk . Let S (respectively S R ) be the homotopy path of H (resp. H R ) which belongs (X sk , 0) (resp. (X sk , 0, k sk )). Denoting p(S R ) the projection of S R on its 10 first coordinates, we have S = p(S R ).
Proof Let f k : R 9 → R be the function defined for all X ∈ R 9 that maps a figure X of K 3,3 to the length between its points P 0 and P 2 and consider the set (iii) ∀(X, t, f k(X)) ∈ S , H R (X, t, f k(X)) = 0, we have S ⊆ S R then S ⊆ p(S R ). Consider now the set p(S R ). Since it is connected, (X sk , 0) ∈ p(S R ) and ∀(X, t) ∈ p(S R ), H(X, t) = 0, we have p(S R ) ⊆ S.

Remark 2
The system of equations (H R ) : H R (X, t, k) = 0 can be rewritten as: . . .
where H (X, t, k) = 0 denotes the system of the 9 first equations involving unknowns (X, t, k).
Considering k as a parameter, (H ) : H (X, t) = 0 is triangular by blocks: each block (B i ) of m i equations involves exactly m i unknowns that do not appear in blocks (B j ) with j < i. Solving one block leads to construct one point.
We will say that (H R ) is quasi-triangularized since except for the d last equations, here d = 1, H R is triangular by blocks. Given G = C[X, A], finding d additional constraints in order to build H R such that (H R ) is quasi-triangularized can be achieved by combinatorial reasoning thanks to the reparameterization method (see [6]). Homotopy function H corresponds to the homotopy function of a new GCSP, say G = C [X, A ], that differs from G only for d constraints. For the example K 3,3 , G is depicted in Fig. 7. Parameters of additional constraints of C \ C (here k) are called driving parameters. Finding the minimum number of additional constraints such that (H R ) is quasi-triangularized is still an open problem.

Exploiting triangularization of (H )
Knowing value for driving parameter k and interpolation parameter t, blocks of (H ) can be solved sequentially to obtain all its solutions. Here each block (B i ) corresponds to a circle-line or a circle-circle intersection, and admits at most two solutions. Since (H ) is usually evaluated in a backtracking process, choosing a solution of (H ) amounts to choosing a branch in a binary tree where level n contains solutions for block (B n ). We can thus index different sequences of choices (i.e. different branches) and note them b 1 , . . . , b n .
For a given branch b ∈ {b 1 , . . . , b n }, we consider two functions: one for the sequential part (H ) of (H R ) and another one for the d last equations. In our example, we define function P b that maps values of t, k to the solution X of (H ) according to branch b. Next, function f : (X, t) → P 2 P 5 − a 8 (t) is for the last equation. Then, solving (H R ) for a given value of t amounts to seek values for Obviously, functions P b are not defined for any values of (t, k) since some intersections may be empty. The domain of definition of P b will be denoted by DP b . It is in general neither open nor closed in the usual topology, but admits a boundary. Here, for instance, points (t, k) such that k + a 1 (t) = a 0 (t) belongs to this boundary.
We will show in section 3 how to use constructive geometry to obtain a quasi-triangularization and continuous functions P b that will be called construction plan. Their domains of definition will be studied and their boundaries characterized in terms of geometric configurations.

Paths of h b
Considering h b as homotopy functions, their homotopy paths can be characterized as one dimensional manifolds with boundaries. We will show in section 4 that the path S R of H R to which belongs the sketch projects itself in a (finite) sequence of paths (S i ) 1≤i≤n for n branches ((S i ) for short), each S i being a path of h bi and having its boundaries in boundaries of DP bi . Hence a point (t, k) belongs to the boundary of S i if and only if figure P bi (t, k) presents a special geometric configuration. Fig. 5 shows the projection on its two last coordinates t and k of a path S R of H R for the GCSP K 3,3 . This closed curve is a sequence of 4 paths S i , which boundaries are represented by black circles.
Considering remark 1, following the homotopy path S of H containing the sketch can be achieved by following this sequence (S i ) of paths, which requires at each prediction step and each Newton-Raphson iteration the inversion of a matrix of size d × d, where d is the number of driving parameters (in the example K 3,3 , d = 1).
However a major difficulty remains: path tracking processes are not robust near boundaries of DP bi . First because h b is not defined anywhere. When approaching a configuration of tangency  (circle-line or circle-circle intersection), the prediction can be made outside DP bi , where there is no intersection. So, no correction is possible. Second, near tangency configurations h b has very high derivatives values, leading quasi certainly the path tracking process to fail.

Changing reparameterization during the path-tracking process
We propose in section 5 a novel approach to overcome this problem that consists in changing the quasi-triangularization H R during the path-tracking process when current point (k, t) is close to the boundary of h b . Suppose the configuration (a) of Fig. 5 is reached when tracking the path S 1 of h b1 from the sketch (sk) = (k sk , 0). It is close to the boundary configuration (b) what may lead the path tracking process to failure.
We change the quasi triangularization H R into H R * by replacing block B 2 by block B 2 * defined as: Parameters: length h0, ..., h7, k Constraints:  where Q is a new point which coordinates are fixed. H R * shares the same properties than H R hence it can be used to define a new construction plan and new homotopy functions h * b . The path S R * of H R * containing the sketch projects itself on its two last coordinates k * and t into a sequence (S * i ) of homotopy paths of h * b . Exploring path S is then reduced to track paths (S * i ) from (a) instead of paths (S i ).
If Q is conveniently positioned (see for instance right part of Fig. 6) current point of tracked path is no more close to a boundary configuration and the tracking process can be continued.
The same method is applied each time a proximity of a boundary configuration is detected. It can lead to introduce new driving parameters when the block B i to be changed does not already involve any driving parameter.

Reparameterized construction plans
One way to get a quasi-triangularization is to perform a geometric solving by a ruler and compass construction. To go on with this geometric point of view, we will use the classical terms of geometric solving. So homotopy paths apply on a construction plan and blocks of equations are instructions.

Construction plan
Locus intersection method (LIM) is a simple method of geometric solving where the entities are sequentially constructed by intersecting loci. These intersections are placed in a list called construction plan. But some problems in 2D and most 3D problems can not be solved by this method. To remedy this, [6] proposes to delete some constraints and add other constraints to make the problem solvable by LIM. A numerical method then seeks to the added parameter values to satisfy the deleted constraints. How to choose the added and removed constraints is non-deterministic and described in [6]. This reparameterization process will not be discussed here.
Consider the GCSP depicted in Fig. 7 that has been obtained from the K 3,3 by substituting distance of parameter h 8 by distance between P 0 and P 2 of parameter k. We will call this GCSP quasi − K 3,3 . This new distance k allows a simple ruler and compass construction for this problem. It consists of a sequence of intersections: (I 0 ) P 0 lying on a line l 0 .
(I 1 ) Construct P 1 as intersection of l 0 and the circle of radius h 0 centered in P 0 .
(I 5 ) Construct P 4 as intersection of circles respectively centered in P 0 and P 4 of radius h 5 and h 4 .
Instruction I 0 fixes a reference to construct figures up to rigid motions. Each other instruction I i constructs a geometric object by intersecting geometric loci defined by both parameters and previously constructed objects. Such a sequence of instructions is called a Construction Plan (CP) of the solutions of quasi − K 3,3 .
More formally, a construction step I i is a term I i [x, A] depending on a set of parameters A and constructing a geometric object x. For instance, the instruction I 2 [P 2 , A 2 ] is written as P 2 = InterCC(P 0 , k 0 , P 1 , h 1 ) where A 2 = {P 0 , k 0 , P 1 , h 1 } is the set of parameters, and P 2 the constructed object.

Definition 3 A Construction Plan P of objects X and parameters
If G = C[X, A] is a GCSP, we say that P = I[X, A] is a construction plan of its solutions if constructed figures fulfill all constraints of C. Fig. 8 shows a CP of solutions of quasi − K 3,3 with reference {P 0 , l 0 }.

Evaluation of a construction plan
For given values A of parameters and A 0 of reference, P = I[X, A] can be evaluated, what means sequentially applying instructions of I. A multi-function, called interpretation, is associated to each instruction of P . For instance, InterCC is linked to a circle-circle intersection procedure. This latter could yield one, two or no solutions (we assume that no solutions are returned if intersection is empty or if circles are the same).
Evaluating a CP leads to an interpretation tree. As for particular equation systems of section 2, each numerical solution is brought by a branch of that tree.
Right part of Fig. 8 shows a sub-tree of the interpretation tree associated to the CP presented in the left part, and the two figures resulting of its evaluation on the two branches in solid line.
Let P = I[X, A] be a construction plan with reference A 0 and I = {I 1 , . . . , I n }.
We note P b (A) = P n b (A) and we call it the numerical function associated to P on branch b. We note P b (A) the figure constructed by this function for values of parameters A.
It is straightforward that if each I i b is smooth (C ∞ ) in its variables A i , P b is smooth in its variables A, as a combination of smooth functions.

Reparameterized construction plan
Let G = C[X, A] be a GCSP and G = C [X, A ] be another one obtained by reparameterization, and P a construction plan for G . C \ C is the set of d added constraints, depending on d driving parameters A + , and C − = C \ C = {c 0 , . . . , c d−1 } is the set of removed constraints.
We call Reparameterized Construction Plan (RCP) the pair R = (P, C − ). If b is a branch of P , we will note P b (A + , A) the numerical function associated to P on branch b to make appear the different roles played by driving parameters and other parameters.
For the following, it is worth mentioning that a RCP is not unique for a given GCSP. By considering constraints of C − with there numerical interpretations, we will note c i (X, A) the value of the numerical interpretation of constraint c i .
We associate to R = (P, C − ) the numerical functions where b is a branch of P , R l is the set of all possible values for A.

Proof of the approach
Assuming that removed constraints can be evaluated for all X, which is the case for distance constraints, solving G for the values A so is equivalent to finding couples (

Border of the domain of definition of a RCP
Be given a branch b, we are interested in characterizing the domain of definition DR b of function R b and more precisely its border ∂DR b , that is a subset of R d × R l . We will show that a point (A + , A) belongs to this border if and only if the figure P b (A + , A) presents one among a finite set of specific geometric configurations.
We recall that the border ∂DR b of DR b is the subset of R d × R l such that for each x ∈ ∂DR b and each ∈ R + , a ball centered in x with ray contains both points of DR b and points of ( Assuming that functions c 0 , . . . , c d−1 are defined for all possibles figures X, DR b is exactly the domain of definition DP b of P b , and ∂DR b is exactly ∂DP b .

Border of DP b
Consider the definition of the numerical function P b given in expression (10). Note DI i b the domain of definition of numerical function associated to instruction Proposition 1 Let P be the sequence of instruction I = {I 1 , . . . , I n } and suppose that each instructions I i involves at least one parameter that is not involved in other instruction. Then following sentences are equivalent.
We prove here that (ii) ⇒ (i) in the case where instruction I i is the circle-circle intersection x i = InterCC(x j , a j , x k , a k ). Suppose a j is a parameter that is not involved in other instructions, Let us enumerate the cases where (x j , a j , x k , a k ) belongs to ∂DI i b ; all these cases correspond to specific geometric configurations.
(1) the two circles are tangent because the triangular equality is satisfied, (2) the two circles are concentric with rays equal to zero, (3) the two circles are concentric with rays not equal to zero.
In cases (1) and (2), we have (A + , A) ∈ DP i b since the two circles have an intersection. A ball centered in (A + , A) of any ray ∈ R + contains values for which the two circles have no intersections: one can simply grow or decrease a j by a sufficient small value. Then (A + , A) ∈ ∂DP i b . In case (3) we consider that I i b is not defined in (x j , a j , x k , a k ), and as a consequence P i b is not defined in (A + , A), since intersections of the two circles are in infinite number. Suppose that one of the centers, say x j , is defined by a previous instruction I j , that involves the free parameter a l . Since I j b is continuous, for any ∈ R + , it is possible to find a point belonging to the ball centered in (A + , A) with ray , resulting from a small variation of a l , such that the two circles intersected in instruction  When I i is a circle-circle intersection, these configurations are listed above. If a point (A + , A) is such that figure P

Border and boundary configurations
When We call boundary of DP b , and we note it ðDP b the set DP b ∩ ∂DP b . We call specific configurations (1), (2) and (4)

Path tracking using RCP
Let G be the GCSP C[X, A], H its associated homotopy function and R = (P, C − ) a RCP involving d driving parameters. This RCP is a quasi-triangularization for the system of equations associated to G. We define the R-extended homotopy function H R : R m × R × R d → R m × R d which is obtained by appending d numerical functions corresponding to added constraints of R to components of H.
Let us recall a R-extended homotopy function for K 3,3 : where a is an interpolation function from A sk to A so , and P 0 P 2 − k corresponds to the added constraint k = distance(P 0 , P 1 ). Note p : R m+1+d → R m+1 (respectively p : R m+1+d → R d+1 ) the projection of an object of R m+1+d on its m + 1 first (resp. d + 1 last) coordinates, and extend these notations to subsets of R m+1+d .
We first generalize the remark 1 of section 2: Proposition 2 If 0 is a regular value of H and numerical functions associated to added constraints are C ∞ over DH × R d then 0 is a regular value of H R . Moreover, if S is an homotopy path of H, it exists a unique homotopy path S R of H R such that p(S R ) = S, and S and S R have the same topology.
This proposition states that following a path S of H can be achieved by following a path S R of H R . The purpose of this section is to show that following S R , and hence S, can be reduced to follow the set p (S R ), which is a one dimensional set in a space of dimension d + 1 instead of Inria Figure 9: A projection of the path S R of the R-extended homotopy function for K 3,3 and its projection p (S R ). Black points mark points of curves that present a boundary configuration. m + 1. However, p (S R ) is in general not a smooth manifold of dimension 1, but we will show that it can be decomposed into a finite sequence of such manifolds that are homotopy paths of functions build from R, and a theoretical algorithm that tracks this sequence will be given. Unfortunately the one who wants to implement this algorithm with a numerical path-tracking method will face difficulties, that are inherent to the multi-functional nature of the numerical interpretation of a construction plan.
From this proposition, two important remarks must be made.
Remark 3 Be given a homotopy path S of H, the path S R in proposition 2 is the subset of R m+1+d whose points are (X, t, A + ) where (X, t) ∈ S and A + are values of driving parameters that can be read on X.
Remark 4 Be given a homotopy path S R of H R and a point (X, t, A + ) ∈ S R , it exists at least one branch b s.t X = P b (A + , a(t)). As a result, be given a point (t, A + ) of p (S R ), it exists at least one branch b s.t X = P b (A + , a(t)) and H R (X, t, A + ) = H(X, t) = 0.
The branch b in remark 4 is unique in most cases. However, if X presents a boundary configuration, there exists at least two such branches in our framework.

R-extended homotopy path
We give here some properties of a connected component S R of (H R ) −1 (0). First, by proving proposition 2, we justify that S R is a one-dimensional smooth manifold. Then we will show that it contains a finite number n of points (X, t, A + ) s.t. X presents a boundary configuration of R.
As a result, one could scan S R as a finite sequence (S R i ) 0≤i<n of pieces of path, which two consecutive elements S R i and S R i+1 share one point (X i , t i , A + i ) s.t. X i presents a boundary configuration of P .
If n = 0, the sequence (S R i ) 0≤i<n is reduced to one element S R 0 = S R . Otherwise, if S R is diffeomorphic to a circle, n = n and each S R i is diffeomorphic to  Fig. 9 presents a projection on its three last coordinates x 9 , t, k of the path S R (that is diffeomorphic to a circle) for the K 3,3 . Black circles mark points that correspond to boundary configuration.

Proof of proposition 2
Let us write the i-th component of H as H i (X) − a i (t), for 0 ≤ i < m. For the sake of the proof, we re-order components of H R : its m first components are those of H and its d last correspond to added constraints and can be written H m+i (X) − k i for 0 ≤ i < d where k 0 , . . . , k d−1 are the driving parameters of R. Hence, noting J H the jacobian matrix of H, the jacobian matrix of H R can be written: is the unique homotopy path of H R s.t. p(S R ) = S and has the same topology than S.

Boundary configurations
Proposition 3 S R contains only a finite number of points (X, t, A + ) s.t. X presents a boundary configuration of R.
When P involves only InterCC and InterCL instructions, boundary configuration are situations (1), (2) and (4) enumerated in subsection 3.3. Given a CP, it exists only finitely, say l, many different boundary configurations. Moreover, each of these can be characterized by a geometric constraint written as a polynomial equation g i = 0.
Consider for instance the instruction P 2 = InterCC(P 0 , k, P 1 , h 1 ) of the RCP for K 3,3 . Configurations (1) or (2) hold for this instruction when circles of respective centers P 0 , P 1 and rays k, h 1 are tangent. It is equivalent to say that P 0 , P 1 and P 2 belongs to the same line, hence it can be translated in the equation det(( − −− → P 0 P 1 , − −− → P 0 P 2 )) = 0. Noting Hg i : R m+1 → R m+1 the function obtained by appending g i to H, a figure (X, t) presents a boundary configuration if and only if there is i s.t. (X, t) is a zero-value of Hg i . We make the following assumption : Next, by lemma on page 11 of [14], we know that roots of Hg i are isolated points. Since Hg i are polynomials, there are finitely many roots.
Cases when this assumption is false corresponds either to cases when parameters and interpolation function are non-generic, or to situations where g i is a consequence of constraints of G, what can be very hard to detect.

Homotopy functions using RCP
Consider now S R as a finite sequence (S R i ) 0≤i<n . Tracking p (S R ) amounts to follow a finite sequence (S i ) 0≤i<n where S i = p (S R i ) for 0 ≤ i < n. Each extremity of S i corresponds to a figure presenting a boundary configuration, and two consecutive elements S i and S i+1 shares such a point.
Even if p (S R ) is not always a manifold and can present for instance self-intersections, we will characterize each S i as a homotopy path of a homotopy function build from R. Fig. 9 presents the path S R for the GCSP K 3,3 , and the projection p (S R ) which each elements S i lie between two black circles.

Definition
We define the homotopy functions where b is a branch of P and R b defined as in equation (11). h b is C ∞ over its domain of definition Dh b = DR b . The border ∂Dh b will be studied. Since functions c 0 , . . . , c d−1 are C ∞ for all (X, t) ∈ R m × R, Dh b (resp. ∂Dh b ) is exactly the domain of definition (resp. the border of the domain of definition) of the function P b • a : (t, A + ) → P b (A + , a(t)).

Boundary configurations
We suppose that the interpolation function a is sufficiently generic and that proposition 1 can be extended to function P b • a. Hence (t, A + ) belongs to ∂Dh b iff figure P b (A + , a(t)) presents a specific geometric configuration. More, supposing that instructions of P are of type InterCC or InterCL, (t, A + ) ∈ ðDh b iff the figure P b (A + , a(t)) presents one of the boundary configurations enumerated in 3.3.
In the latter case, it exists b = b s.t. (t, A + ) belongs to ðDh b and P b (A + , a(t)) = P b (A + , a(t)), and finally h b (t, A + ) = h b (t, A + ). Moreover, we suppose that b is unique. It is not when X presents a boundary configuration for more than one instruction, but assumption 1 rules these cases out.
In other words, if (t, A + ) ∈ ðDh b , then (P b (A + , a(t)), t, A + ) is a solution of H R = 0 with multiplicity two, and two branches b and b allow to construct P b (A + , a(t)).

h b vanishes on p (S R )
Consider a point (t, A + ) of p (S R ). It exists a branch b of P s.t. the point (P b (A + , a(t)), t, A + ) belongs to S R . Hence h b (t, A + ) = 0. More, one can state: Let 0 ≤ i < n. S i is connected since it is the projection of S R i . Consider two points (X 1 , t 1 , A + 1 ) and (X 2 , t 2 , A + 2 ) of S R i and suppose that X 1 does not present a boundary configuration. Hence it exists a unique branch b s.t. X 1 = P b (A + 1 , a(t 1 )) and h b (t 1 , A + 1 ) = 0. Since the segment of path S R i between (X 1 , t 1 , A + 1 ) and (X 2 , t 2 , A + 2 ) is a smooth deformation of (X 1 , t 1 , A + 1 ) into (X 2 , t 2 , A + 2 ) that does not pass trough a boundary configuration, we have X 2 = P b (A + 2 , a(t 2 )) and h b (t 2 , A + 2 ) = 0.

Homotopy paths of h b
We note h b|∂Dh b the restriction of h b to the boundary of Dh b , and we adapt to our notations a lemma that can be found on page 13 of [14]: Assuming this regularity condition, components of h −1 b (0) are either diffeomorphic to circles or lines that do not have any boundary, or to segments [0, 1[ or [0, 1] that have boundaries. In these last cases, a boundary point (t, A + ) is such that P b (A + , a(t)) presents a boundary configuration, and it exists a unique b = b with: Hence, considering proposition 4, for each sub-set S i of p (S R ), there is a branch b i s.t. S i is a homotopy path of h bi with extremities in ðDh bi and can be followed by a path-tracking method since it is a one dimensional manifold.

Regularity of h b
It seems difficult to transfer properties of regularity of value 0 from H and Hg i , for 0 ≤ i < l, to functions h b and h b|∂Dh b , for all branches of P . Instead of showing that they hold, and then that hypothesis of lemma 1 are satisfied, we simply give some tracks to ensure that sets S i are sufficiently regular.
First, S i does not present any self-intersection. Consider S R i and the branch b such that X = P b (A + , a(t)) for (X, t, A + ) ∈ S R i . If (t, A + ) is a node of S i then (X, t, A + ) is a node of S R i what contradicts the result of proposition 2. Remark that p (S) can present self-intersections, but in this case it is an intersection of sets S i and S j with i = j. Then, S i \ ðS i is smooth: instead of a projection, it can be seen as the image of pieces of S by the smooth function that maps to a figure values of driving parameters.
Finally, we assume that cases where the tangent of S i \ðS i vanishes corresponds to non-generic cases.

Path tracking
Let us now expose the main result of this section. It states that an homotopy path S of H can be deduced from a sequence of homotopy paths of homotopy functions h b , build from a RCP.
We recall that following S requires to invert the jacobian matrix of H which size is m × m, while following homotopy paths of h b leads to invert matrix of size d × d where d m. The possible gain in computation time brought by this approach justifies our motivation.

Proof of proposition 5
Let S be an homotopy path of H, and 0 be a regular value of H. From proposition 2 it exists a unique path S R satisfying p(S R ) = S and from proposition 3 the latter can be parsed in a finite sequence (S R i ) 0≤i<n . The projection p (S R ) is then parsed in a finite sequence (S i ) 0≤i<n which each element S i is an homotopy path of a function h bi . More, S i has its boundaries in ðDh bi and is a one dimensional manifold that can be followed by a path tracking method, thanks to lemma 1. The sequence ((S i , b i )) 0≤i<n clearly satisfies points (i),(ii) and (iii).

Theoretical path-tracking algorithm
The above proposition gives rise to a theoretical path-tracking algorithm that allows to explore a homotopy path S ∈ R m+1 of H by following a sequence of path (S i ) 0≤i<n with S i ∈ R d+1 . Continuation solving of RCP Input: GCSP C[X, A], values A so , sketch X sk , RCP R = (P, C − ) with driving parameters A + . Output: list of solutions L sol 1. Read values A sk and A + sk from sketch X sk , build an interpolation function a from A sk to A so . a(t)). • X = P b (A + , a(t)) converges to a non-derivable figure, then stop;

Find the current branch
• (t, A + ) crosses hyperplane t = 1, then append solution X = P b (A + , a(t)) to L sol ; • (t, A + ) ∈ ðDh b , then find the unique branch b s.t. P b (A + , a(t)) = P b (A + , a(t)), do b = b and go to 4.

Implementation with a numerical method
Regardless to the fact that it fails if assumption 1 is not satisfied, the one who wants to implement the above algorithm with a numerical path tracking method will face two major difficulties. The first one lies in the nature of numerical methods: calculated points are approximations of paths S i , and (t, A + ) ∈ ðDh bi will never be exactly crossed. Hence a prediction-correction method will terminate with a predicted point outside of Dh bi for which it is not possible to apply the correction method. A piecewise linear approximation will terminate when computing a simplex that has a vertex outside of Dh bi . Linear approximation of h bi can not be computed in this case.
This issue could be addressed by the use of commonly called watch-dogs to estimate to which geometric configurations correspond figures that belong to the set {P bi (A + , a(t))|(t, A + ) ∈ ðDh bi } in order to identify the branch b i+1 such that ðDh bi = ðDh bi+1 . Then, by appending an equation describing this geometric configuration to h bi , a point approximating S i ∩ ðDh bi can be obtained by applying a numerical method, for instance Newton-Raphson iterations, to the obtained numerical function.
The second difficulty comes from the fact that more the current approximation (t, A + ) of S i is close to the boundary of Dh bi , higher are the values of the derivatives of h bi .
It can be easily remarked by considering a function that maps to a positive real number r the positive y-coordinate of the intersections of the unit circle with a circle centered in (0, 2) of radius equal to r. Noting y this function, we have y(r) = √ −r 4 + 10r 2 − 9 4 and y (r) = −4r 3 + 20r 8 √ −r 4 + 10r 2 − 9 . It is straightforward that Dy = [1,3] and that ðDy = {1, 3} corresponds to configurations where the circles are tangent; we have lim r→∂Dy y (r) = +∞. Notice that it is not due to the chosen system of coordinates. It can gives rise either to a growing inexactitude when approaching a boundary configuration if the path is tracked with a constant step. Otherwise, if the step is adapted during the tracking, the former can be decreased infinitely.

On-the-fly change of RCP
We propose here a new approach to explore entirely a homotopy path S of H constructed from a GCSP G without tracking it, that uses reparameterization to decrease the cost of path tracking, but overcomes the difficulties pointed out in section 4.
As stated before, a reparameterization is not unique. Moreover, when a figure presents a boundary configuration for a given reparameterization, the latter can be slightly modified in order to avoid this situation. Subsection 5.1 shows how to find this new reparameterization. An algorithm will be proposed in subsection 5.2 which basic idea is to change on-the-fly, during the path tracking process, the reparameterization used to define the tracked path when it is detected to be close to a border configuration. The implementation of this algorithm requires to specify the notion of closeness to a boundary configuration. It is done in subsection 5.3 by defining a dedicated geometric criteria.
In this section, G = C[X, A] is a 2D GCSP. For the sake of simplicity it involves only constraints of distances. H is its associated homotopy function.

Different RCP
For a given problem G, we put forward two points about considering different RCP which our method will take advantage. Let S be a homotopy path of H.
The following assertion that comes directly from proposition 2 and states that different RCP lead to follow the same path.
Proposition 6 Let R 1 and R 2 be two RCP of G, then it exists a unique homotopy path S R1 (resp. S R2 ) of the R 1 -extended (resp. R 2 -extended) homotopy function such that p(S R1 ) = p(S R2 ) = S.
As a result, applying the procedure Continuation solving of RCP described in section 4 with different RCP of G leads to compute projections of the same path S.
But it is important to notice that different RCP have different boundary configurations.
We prove this proposition by showing how to construct R 2 . Let R 1 = (P, C − ), with d driving parameters A + , and P = I[X \A 0 , A∪A 0 ∪A + ] (A 0 are the unknowns used to build the reference). Since G involves only distance constraints, instructions of P have type InterCC or InterCL.
If X presents a boundary configuration, it exists a smaller k s.t. X presents a boundary configuration of instruction I k of P . Consider first that I k does not involve any driving parameter.
Let I k be x k = InterCC(x i , a i (t), x j , a j (t)); hence in the figure X, circles of centers x i and x j with radii a i (t) and a j (t) are either tangent (case a), or concentric with vanishing radii (case b).
In case a, consider a new point y = x i , x j . Then circles of centers x i and y, with radii a i (t) and ||yx k || 2 are not tangent, and intersect in x k , as it is illustrated in the left part of Fig. 10.
One can define the instruction I k as x k = InterCC(x i , a i (t), y, b) and the construction plan for which X does not present a boundary configuration.
In case b, it is necessary to introduce two new points y and z and two new driving parameters to construct the instruction I k .
The construction of I k when I k is an InterCL instruction is illustrated in right part of Fig.  10.
Consider now that I k involves a driving parameter r, say x k = InterCC(x i , a i (t), x j , r), and that circles of centers x i and x j with radii a i (t) and r are tangent (the other cases are similar). y being a new point defined as in case a, we define I k as x k = InterCC(x i , a i (t), y, b) with b a new driving parameter, P as previously, and R 2 as [P , C − ]. Hence R 2 is a RCP of G with d driving parameters (A + \ {r}) ∪ {b} for which X does not present a boundary configuration.

Changing RCP during the path-tracking process
We rely on propositions 6 and 7 to design a new algorithm to compute the path S of H containing the sketch using reparameterization to decrease the cost of the tracking process. To each instruction I k of the construction plan is associated a geometric criterion γ k which role is to decide when a figure X corresponding to a point of the tracked path is close to a boundary configuration. When this situation arises, the RCP used to define the tracked path is modified as shown in the proof of proposition 7. This new RCP defines a new path that is followed from the point corresponding to the last computed figure, until it is close to a new boundary configuration (and the RCP is one more time changed) or it reaches a termination condition.
Applying this method leads to track a sequence of pairs (S i , R i ) i∈N ⊆N where N is an internal of N and for each i, S i is a homotopy path of the homotopy function using RCP R i . In 5.2.3 we will show that if the algorithm terminates, i.e. N is bounded, the set of computed figures is S. Subsection 5.2.4 justifies termination.

Main algorithm
To each instruction I k , for 1 ≤ k ≤ n, we associate the real univariate C ∞ function γ k , that maps to a figure X an element of [0, 1]. Be given α ∈]0, 1[, we consider that a figure X is close to a boundary configuration when min k γ k (X) ≤ α. Functions γ k are defined in subsection 5.3.
Continuation solving of RCP with adaptive CP Input: G = C[X, A], values A so , sketch X sk , RCP R = (P, C − ) with driving parameters A + , α ∈]0, 1[. Output: list of solutions L sol 1. Read values A sk from sketch X sk , build an interpolation function a from A sk to A so ; let t = 0.
3. Perform procedure Initialize RCP with R and X as input values.

Read values A +
cur on X cur , find the current branch b s.t. P b (a(t), A + cur ) = X cur and build homotopy function h b (A + , t) = R b (A + , a(t)). (t), A + ) loops on X sk while checking if :

Follow path of
• X converges to a solution at infinity, then stop; • X converges to a non-derivable figure, then stop; • (A + , t) crosses P 1 , then append solution X to L sol .
• It exists I k [X k , A k ] in P s.t. γ k (X) < α, then perform procedure Change RCP with R, X, k and α as input values, let X cur = X and go to 4.
Steps 1 to 3 are initialization steps. The procedure Initialize RCP is rather technical, and is informally described below; its necessity will emerge from the procedure Change RCP.
Step 4 is performed each time the RCP is changed. It consists in reading on X cur values of new driving parameters A + cur and branch b. Suppose the instruction introduced when changing the RCP is x k = InterCC(x i , a i (t), y, s), then s is set to the value of the distance between x k and y. Finding b is simply performed by evaluating P and choosing at each step the result among that two corresponds to X cur . The cost of this step is linear in the number of instructions of P .
A technical point lies in step 5: how to ensure that there is no track-back on S between S i and S i+1 . It is addressed by checking the progression on S.

Procedure Change RCP
To change the RCP in the procedure Change RCP we consider the instruction I k which current figure X reaches the boundary of the domain of definition. That is we have min i γ i (X) = γ k (X) < α.
This procedure introduces new driving parameters, and we need to distinguish among the set of driving parameter those that were involved in the RCP given as input of the main procedure. They will be called original.
Suppose I k is x k = InterCC(x i , a i , x j , a j ), and switch on the number of driving parameters in {a i , a j }.

Inria
Case 0 One couple among (a i , x i ), (a j , x j ), say (a j , x j ), is chosen to construct I k and the con- where y is a new point which placement is described in 5.3 and b a new parameter. R is modified as follows: , and the procedure is stopped.
Case 1 Let a j be the driving parameter: either it has been introduced at the time of a previous call by case 0, then perform step a), or it is an original driving parameter then perform step b).
a) Let I k be the instruction x k = interCC(x i , a i , x j , a j ) that was replaced by I k .
We assume x j ∈ A 0 . Since it is a reference, it can be moved in combination to a j without modifying X. Then x j is moved in such a way that γ k (X) ≥ α and a j is read on X ∪ {x j }.

Case 2 This situation never happen by construction.
This procedure raises three observations. First it is assumed in step 1b) that if a j is a driving parameter in x k = InterCC(x i , a i , x j , a j ) then x j ∈ A 0 , what is not true for original driving parameters. Actually, the role of the procedure Initialize RCP is to introduce a new point y, with y = x j , for each driving parameter of the original RCP, and to appropriately modify the instruction involving a j .
The second point is that procedure Change RCP increases the number of driving parameters each time it enters in case 0, and decreases it only if it performs step 1a) when the original instruction can be restored. To decrease more systematically the number of driving parameter, an initial step is perform when entering procedure Change RCP. It consists in applying Remove driving parameters with R, X, k and α as inputs. The procedure Remove driving parameters accepts as inputs a RCP , the current figure X and the index k of the instruction to be modified. For each instruction I j with j = k involving a non-original driving parameter, it restores the corresponding original instruction I j if γ j (X) > α.
The third observation is that it is straightforward to transpose this procedure to treat interCL instructions.

Correctness
The goal of this section is to show that when applying procedure Continuation solving of RCP with adaptive CP, the projection in the figures space of the obtained sequence of paths is S, the path of H containing the sketch X sk . Let (S i , R i ) i∈N ⊆N be the sequence of couples (paths,RCP) computed when applying procedure Continuation solving of RCP with adaptive CP, and suppose N is bounded. projection of a point (or a subset) of R m × R × R di on its d i + 1 last coordinates. Recall that p is the projection on the m + 1 first coordinates.
Each path S i is the projection by p i of a 1-dimensional manifold S Ri that is a subset of a homotopy path of H Ri . Consider (S Ri ) i∈N . We will show that if S is diffeomorphic to a circle, i∈N p(S R1 ) = S.
When performing for the first time the step 4 one have X cur = X sk , hence noting P the construction plan of R 0 , b is s.t. P b (A + sk , 0) = X sk , then (X sk , 0, A + sk ) ∈ S R0 and we deduce that p(S R0 ) ⊆ S.
If N = {0}, i.e. procedure terminates without entering in Change RCP because path overlaps, then S R0 is diffeomorphic to a circle and p(S R0 ) = S.
Otherwise, naming X cur the figure corresponding to the last point (t i , A i ) of S i computed before entering the procedure change RCP for the i + 1-th time, P the construction plan of R i+1 and A + its driving parameters, then when performing for the i + 1-th time the step 4, b is chosen s.t. P b (A + sk , t i ) = X cur . Then p(S Ri ) and p(S Ri+1 ) are connected subsets of a homotopy path of H. If there is no backtrack, p(S Ri ), p(S Ri+1 ) share exactly one point if N contains more than two elements 3 . Considering p(S Ri ) ⊆ S as a recurrence hypothesis, we have p(S Ri+1 ) ⊆ S.
If the figure corresponding to the last extremity of the tracked path is the sketch, then i∈N p(S Ri ) is diffeomorphic to a circle, what finishes the proof. When S is diffeomorphic to an interval the procedure has to be performed two times, to explore S from the sketch in both directions.

Termination
Here we consider a path tracker as an abstract object that can be iterated. It either successes and progresses along a path, returning a new point and stops when a termination condition is reached or fails only if current figure X reaches a boundary configuration, i.e. min i γ i (X) = γ k (X) ≤ α.
Applying our algorithm yields to perform a sequence s of steps, each being of one of the types: In a step of type (a), a termination condition is eventually reached and the algorithm stops without error.
We build our discussion on a disjunction on infinite sequences s knowing that a step (a) is followed either by a step (a) or a step (b), a step (b) is followed by a step (c), and a step (c) is followed either by a step (a) or (b).
Suppose s contains an infinite number of consecutive (b) and (c) steps: the procedure change RCP is inefficient. A straightforward necessary condition that avoid this situation is the following. If X is the last figure computed, R i the RCP which X satisfies γ k (X) ≤ α and R i+1 the new RCP constructed by change RCP then γ k (X) > α if γ k is the geometric criterion associated to the new instruction I k in R i+1 . In subsection 5.3 are given a definition of γ k and a geometric placement of new points that ensure this necessary condition. Assuming this hypothesis, the RCP is invariant after a finite number of consecutive steps (b) and (c), and the current figure X satisfies min i γ i (X) > α. Hence the path-tracker should progress along the path.
Consider s contains an infinite number of consecutive steps (a). Recall that the tracked path S i is the projection of the path S Ri that is diffeomorphic to a circle and has a finite arc-length. We have made the assumption that each iteration of the path-tracker progresses along the path. Measuring this progression with the arc-length of S i between two consecutive points, we suppose that it exists δ > 0 s.t. at each iteration, this progression is greater that δ. Hence a termination condition is reached after a finite number of steps (a).
Suppose s is an infinite sequence of patterns a. . .abc. . .bc, where each pattern is composed by a finite number of steps (a) followed by a finite number of steps (b) and (c). Each iteration of the pattern yields a progression along a path S i , hence it yields a progression along S. Nevertheless, it seams not clear that it exists δ > 0 s.t. each progression along S is greater than δ. Think for instance to a sequence of (S i ) i which arc-lengths grow and are unbounded with i. Remark that for a given GCSP, it exists only a finite number of RCP, and a finite number of R-extended homotopy functions, each admitting a unique homotopy path S R that contains the sketch. The geometry of S R is linked to the placement of points added by the procedure Change RCP. But if we consider that figures X and added points stays in a compact set, S R belongs to a compact and has a bounded arc-length. Finally, considering that it exists a longest path S R , we claim that is exists δ > 0 that is smaller than the progression along S at each iteration. Hence a termination condition is reached after performing a finite number of patterns a. . .abc. . .bc.

Geometric criteria
The main procedure given in 5.2.1 and procedure change RCP need geometric criteria to be implemented. First, for each instruction I k , with respect to its symbol, we associate in 5.3.1 a real positive function γ k that maps to a figure an estimation of how close is the figure to a boundary configuration. Then change RCP, when it performs steps 1. (resp. 2b) introduces a new point (resp. moves a reference point). A way of placing this new point will be given in 5.3.2.
Finally, we will show in 5.3.3 that the placement of the new point is coherent with the stopping criterion, i.e if I k is such that γ k (X) ≤ α and I k is the new instruction then γ k (X) > α.

Stopping criteria
where c is the distance between x k and the line x i x j . If I k is the instruction x k = InterCL(x i , a i , x j ), we define where c is the distance between x k and the projection p of x k on line x j .

Placement of new points
In steps 1 and 2.b of procedure change RCP, one needs to place (i.e. to find coordinates of) the point y.
If I k is the instruction x k = InterCC(x i , a i , x j , a j ), let p be the projection of x k on the line x i x j , and − − → x k p/x k p the unit vector which direction is perpendicular to the line x i x j . We use the placement of y: If I k is the instruction x k = InterCL(x i , a i , x j ), we propose to place y on the line x j , at a distance a i of x k : where p is the projection of x i on line x j .

Coherence of criteria
Let α ≤ 1 and suppose γ k (X) is such that γ k (X) < α. Suppose that I k is the instruction obtained after applying the procedure change RCP, with placement of y as described above. If α ≤ 1/2, then γ k (X) > α.
It is straightforward when I k is an InterCL instruction; see right part of Fig. 11. When I k is c 2a i and to deduce the result. It is illustrated by left part of Fig. 11.

Implementation and results
The procedure Continuation solving of RCP with adaptive CP described in previous section has been implemented in C++, giving rise to a program that accepts a GCSP G, an initial RCP R of G, a sketch X sk and values A so . Numerical path-tracking process is achieved by a simple prediction-correction method with fixed step. It is described in subsection 6.2.
In subsection 6.3, we present numerical results about the application of our method on four point-point distances problems, one in 2D and three in 3D. It confirms that the approach of [9] that consists in following the path S of H which belongs the sketch is fast and provides several solutions (sometimes all of them) of problems that resist to divide and conquer methods and are two large to be solved by classical homotopy or branch and bound approach. The method presented here uses parameterized construction plans and brings an important speed-up of this approach while staying robust.
We first give details about its generalization to 3D cases.

Generalization to 3D GCSP
3D GCSP with point-point distance constraints fit well to the solving method depicted in this paper: results of [9] as well as reparameterization approach stay valid. One important issue when solving a 3D GCSP is to decide if it is generically well constrained since there is not known combinatorial or structural characterization of the latter property. It is assumed to hold for GCSP considered here. (I 2 ) P 2 = InterSSP (P 0 , h 1 , P 1 , h 2 , pl 0 ) where InterSL holds for a sphere-line intersection, InterSSP for the intersection of two spheres and one plane, and InterSSS for a three spheres intersection.
It is simpler to see a three spheres (resp. two spheres and one plane) intersection as the intersection of two spheres, which is a circle of non strictly-positive ray, and then the intersection of a circle and a sphere (resp. a plane). We will write InterSS for two spheres intersections, and InterSC (resp. InterP C) for sphere-circle (resp. plane-circle) intersections.

Border and boundary configurations
Possible border and boundary configurations of 3D parameterized construction plans can be enumerated and are similar to those encountered in a 2D context.
InterSS and InterSL are straightforward generalizations of InterCC and InterCL instructions.
InterSC instructions can be simplified in InterCC instructions by considering the intersection of the sphere and the plane which belongs the circle. InterP C is seen as a InterCL instruction provided that the plane is not parallel or equal to the plane which belongs the circle.

Stopping criteria and placement of new points
To generalize the algorithm in a 3D universe, one need to define stopping criteria and placement of new points for 3D instructions.
For InterSS and InterSL, instructions, we use the direct 3D generalization of stopping criteria and placement of new points of InterCC and InterCL instructions. Consider for instance the instruction x k = InterSS(x i , a i , x j , a j ). Then the stopping criterion γ k (X) associated to it is c/max(a i , a j ) where c is the distance of x k to the line x i x j . Suppose a j < a i . When swapping this instruction with x k = InterSS(x i , a i , y, b), one can choose b = a i and y = where p is the projection of x k on line x i x j . We propose the following definition for the stopping criterion associated to an InterSC instruction. Let x k = InterSC(x s , a s , x c , a c ), x be the projection of x s on the plane pl which belongs the circle (x c , a c ), and c be the distance from point x k to line x c x. Then γ k (X) = c/max(a s , a c ).
When entering the procedure change RCP, if γ k (X) < α, the instruction x k = InterSC(x s , a s , x c , a c ) where p is the projection of x k on line x c x and b = a c . It is easy to show that this criterion is coherent in the sense of 5.3.3. The placement of y is depicted on Fig. 12. 6.2 Implementation

Path tracking
A path S is followed from an initial point (t 0 , d + 0 ) by a classical prediction-correction method: prediction is performed by an Euler predictor with constant step δ ∈ R and correction by Newton-Raphson iterations. Notice that it is a naive path tracking algorithm that does not ensure that there is no jump between two close paths. Clever approaches which merely rely on interval arithmetic (see [10,4,13]) have to be used to guarantee a path tracking algorithm. Here we just suppose that δ is small enough to follow considered curve. In the examples considered, such jumps never occurred.
The jacobian matrix J(t 0 , d + 0 ) is numerically computed with the finite difference method, and the 1-dimensional kernel T of the linear application J(t 0 , d + 0 ) is obtained with basics linear algebra tools. A prediction (t 0 , d + 0 ) is then computed in the direction of T , with prediction step δ. This prediction is then corrected with Newton's method in (t * 0 , d + 0 * ) ∈ S such that (t * 0 , d + 0 * ) ∈ T p , where T p is the hyperplane perpendicular to T . When (t 0 , d + 0 ) and the prediction (t 0 , d + 0 ) are such that (t 0 − 1).(t 0 − 1) < 0 (i.e. hyperplane P 1 lies between these two points), a new prediction (t 0 , d + 0 ) such that t 0 = 1, is performed, and (t 0 , d + 0 ) is corrected in (1, d + 0 * ) to precise the solution.
Prediction is made along line T , and the unit tangent vector is obtained by choosing an orientation on T . An orientation is arbitrarily chosen at the first iteration while for next iterations orientations are those that minimize the angle between two consecutive tangent vectors. If the followed path is diffeomorphic to a line, the tracking algorithm has to be applied a second time with the opposite orientation.
When the path-tracking process stops because stopping criterion is satisfied, some driving parameters are added or removed. As a consequence, the next part of the path will not be tracked in the same space. To reinitialize the tracking algorithm, one of two orientations on the line T has to be chosen again. Last used tangent vector is translated in the new space in order to compare it with unit vectors of T .
A loop can be detected when a point of the curve is close enough to the initial point (0, d + sk ). A loop detection test is then performed when hyperplane P 0 lies between previous and current points. In such a case, we search a point on P 0 by using the same mechanism as described above to precise a solution on P 1 and compare it with initial point in the space of original driving parameters.

Differentiation of the CP
When applying the above path-tracking method to track the homotopy path of a function h b , most of the computation time is spent in the evaluation of the underlying CP. Most evaluations intervene in computation of jacobian matrices by finite difference method that needs about d evaluations, and such matrices are computed at each prediction step and at each iteration of the Newton-Raphson method in a correction step.
It makes the differentiation one of the operations to improve in priority. Here we exploit the acyclic computation scheme of a construction plan to improve it. Suppose driving parameters k = (k 0 , . . . , k d−1 ) appear in steps I i0 , . . . , I i d−1 with i d−1 ≥ . . . ≥ i 0 . The vector ∂h b ∂k d−1 (t, k) of the jacobian matrix is obtained by evaluating h b in (t, k 0 , . . . , k d−1 − ) and in (t, k 0 , . . . , k d−1 + ) for a sufficiently small . Remark that figures resulting of these evaluations differs only from the step i d−1 since I i d−1 is the first step that involves k d−1 .

Inria
As a result, a manner of optimizing the cost of the differentiation in (t, k) of h b is to evaluate it entirely a first time in (t, k), to save the obtained figure and hence to compute partial derivatives ∂h b ∂ki in order of decreasing index by evaluating the CP from step I ki . Our implementation incorporates this optimization.

Results
We give here numerical results about application of the path tracking mechanism presented in this paper to point-distance problems, one in 2D and three in 3D. Table 1 presents these results. For each problem, are given the number of equations of the numerical system F (row "Dof/Dor") and running times and number of solutions when applying a classical homotopy method to solve F . It is implemented by the free software HOM4PS-2.0 [18]. Cost in computation time increases dramatically with the degree of freedom of the system due to growing number of path to be followed.
Then row "first path of H" gives the number of solutions found and times needed to follow the homotopy path when homotopy function used is H defined as a system of equations (so with no triangularization).
Finally, the row "Presented method" gives running times of Continuation solving of RCP with adaptive CP. Details about tuning main parameters of this process are presented below before a brief description of its execution for each problem.

Details about parameters
The interpolation function a is defined as follows: a has a compact positive support and fulfills assumptions of characterization of homotopy paths presented in [9]. For each problem, the prediction step δ of the path tracker has been chosen as the biggest such that the process terminates with success. Stopping criterion α has been set to value 0.1 that seems to fit well to our algorithm.

The problems
The goal of the octahedron problem is to construct a solid with 6 vertices, 12 edges and 8 triangular faces knowing the lengths of its 12 edges. A quoted sketch of this problem is given in left part of Fig. 13. This problems is related to the parallel robot called Stewart platform and results in a system F of 12 equations involving 12 unknowns. With parameters: it admits 4 solutions up to rigid motions. These solutions are found in about 27 seconds by HOM4PS2. When using the sketch of Fig. 13, these 4 solutions lies on the path of the sketch, hence they are all found when following it. The second problem comes from molecular chemistry and is picked up from [16]. Coordinates of 8 points in the 3D space have to be found knowing 18 pairwise distances. It corresponds to a disulfide molecule (see right-hand part of Fig. 13). A valuation of parameters is exhibited in [16] that leads to 18 solutions all found by a branch-and-prune algorithm in more than 10 minutes. With these parameters, only 8 solutions lies on the path of the sketch but there are found in 0.4s.
Dodecagon and Icosahedron problems are illustrated on Fig. 14. The former gives rise to a numerical function F with 21 constraints, for which 2 solutions were found both by HOM4PS2 and our approach when using parameters: The numerical system of equations associated to the icosahedron problem involves 30 equations and is too large to be solved by a classical homotopy method. With our approach, we find 28 solutions within 4 seconds. Table 2 presents details about execution of Continuation solving of RCP with adaptive CP for each problem. The rows "Path of H" recall times of tracking paths of function H and give  the number of iterations of prediction-correction needed. The rows "Presented method" recall running times, give the number of iterations of prediction-correction, the number of changes of the RCP (row "RCP changing") and finally the average (and maximum) number of driving parameter of the RCP (average is weighted by the number of prediction-correction performed with), in row "DP: average (max.)". In both cases, row "t/i in ms" gives the average time spent in one prediction-correction iteration. It shows how using parameterized construction plans does speed-up the path tracking process.

Conclusion
Well-constrained GCSP often have many solutions. The existing solvers offer either a single solution or all the solutions. These latter are of limited practical interest because either the class of problems they solve are reduced or their complexity is exponential. But even if not all solutions are needed, an effective method must provide to the user several solutions close to the sketch.
In this article we present a method for solving geometric constraints problems by homotopy that is effective and robust. Regarding efficiency, the method does not follow all paths but only one starting from the sketch. We then find the nearest solution of the sketch but also other solutions. Indeed we are in the conditions of having a homotopy path diffeomorphic to a circle to which belong usually more than one solution.