PTOPO: Computing the Geometry and the Topology of Parametric Curves
PPTOPO : Computing the Geometry and the Topology of ParametricCurves
Christina Katsamaki, Fabrice Rouillier, Elias Tsigaridas
Inria Paris, IMJ-PRG, Sorbonne Universit´e and Paris Universit´e
Zafeirakis Zafeirakopoulos
Institute of Information Technologies, Gebze Technical University, Turkey
Abstract
We consider the problem of computing the topology and describing the geometry of a parametriccurve in R n . We present an algorithm, PTOPO , that constructs an abstract graph that is isotopic tothe curve in the embedding space. Our method exploits the benefits of the parametric represen-tation and does not resort to implicitization.Most importantly, we perform all computations in the parameter space and not in the implicitspace. When the parametrization involves polynomials of degree at most d and maximum bitsizeof coe ffi cients τ , then the worst case bit complexity of PTOPO is (cid:101) O B ( nd + nd τ + d ( n + n τ ) + d ( n τ + n ) + n d τ ). This bound matches the current record bound (cid:101) O B ( d + d τ ) for the problemof computing the topology of a plane algebraic curve given in implicit form. For plane and spacecurves, if N = max { d , τ } , the complexity of PTOPO becomes (cid:101) O B ( N ), which improves the state-of-the-art result, due to Alc´azar and D´ıaz-Toca [CAGD’10], by a factor of N . In the same timecomplexity, we obtain a graph whose straight-line embedding is isotopic to the curve. However,visualizing the curve on top of the abstract graph construction, increases the bound to (cid:101) O B ( N ).For curves of general dimension, we can also distinguish between ordinary and non-ordinaryreal singularities and determine their multiplicities in the same expected complexity of PTOPO byemploying the algorithm of Blasco and P´erez-D´ıaz [CAGD’19]. We have implemented
PTOPO in maple for the case of plane and space curves. Our experiments illustrate its practical nature. Keywords:
Parametric curve, topology, bit complexity, polynomial systems
1. Introduction
Parametric curves constitute a classical and important topic in computational algebra andgeometry (Sendra and Winkler, 1999) that constantly receives attention, e.g., (Sederberg, 1986;Cox et al., 2011; Bus´e et al., 2019; Sendra et al., 2008). The motivation behind the continuous
Email addresses: [email protected] (Christina Katsamaki),
[email protected] (Fabrice Rouillier), [email protected] (Elias Tsigaridas), [email protected] (ZafeirakisZafeirakopoulos) a r X i v : . [ c s . S C ] J a n nterest in e ffi cient algorithms for computing with parametric curves emanates, among othersreasons, by the omnipresence of parametric representations in computer modeling and computeraided geometric design, e.g., (Farouki et al., 2010).We focus on computing the topology of a real parametric curve, that is, the computation of anabstract graph that is isotopic (Boissonnat and Teillaud, 2006, p. 184) to the curve in the embed-ding space. We design a complete algorithm, PTOPO , that applies directly to parametric curves ofany dimension. We consider di ff erent characteristics of the parametrization, like properness andnormality, before computing the singularities and other interesting points on the curve. Thesepoints are necessary for representing the geometry of the curve, as well as for producing a certi-fied visualization of plane and space curves. Previous work.
A common strategy when dealing with parametric curves is implicitization.There has been great research e ff orts, e.g., Sederberg and Chen (1995); Bus´e et al. (2019) and ref-erences therein, in designing algorithms to compute the implicit equations describing the curve.However, it is also important to manipulate parametric curves directly, without converting themto implicit form.The study of the topology of a real parametric curve is a topic that has not received muchattention in the literature, in contrast to its implicit counterpart (Diatta et al., 2018; Kobel andSagralo ff , 2014). It requires special treatment, since for instance it is not always easy to choosea parameter interval such that when we plot the curve over it, we include all the important topo-logical features (Alc´azar and D´ıaz-Toca, 2010). Moreover, while visualizing the curve usingsymbolic computational tools, the problem of missing points and branches may arise (Recio,2007; Sendra, 2002). Alc´azar and D´ıaz-Toca (2010) study the topology of real parametric curveswithout implicitizing. They work directly with the parametrization and address both plane andspace real rational curves. Our algorithm to compute the topology is to be juxtaposed to theirwork. We also refer to Caravantes et al. (2014) and Alberti et al. (2008) for other approachesbased on computations by values and subdivision, respectively.To compute the topology of a curve it is essential to detect its singularities. This is an im-portant and well studied problem (Alc´azar and D´ıaz-Toca, 2010; Rubio et al., 2009; Kobel andSagralo ff , 2014) of independent interest. To identify the singularities, we can first computethe implicit representation and then apply classical approaches (Walker, 1978; Fulton, 1969).Di ff erently, we can compute the singularities using directly the parametrization. For instance,necessary and su ffi cient conditions to identify cusps and inflection points are expressed in theform of determinants, e.g., (Li and Cripps, 1997; Manocha and Canny, 1992).On computing the singularities of a parametric curve, a line of work related to our approach,does so by means of a univariate resultant (Abhyankar and Bajaj, 1989; van den Essen and YU,1997; Park, 2002; P´erez-D´ıaz, 2007; Gutierrez et al., 2002a). We can use the Taylor resultant(Abhyankar and Bajaj, 1989) and D -resultant (van den Essen and YU, 1997) of two polynomialsin K [ t ], to find singularities of plane curves parametrized by polynomials, where K is a field ofcharacteristic zero in the first case and arbitrary in the latter, without resorting to the implicitform. Park (2002) extends previous results to curves parametrized by polynomials in a ffi ne n -space. The generalization of the D -resultant for a pair of rational functions and its application tothe study of rational plane curves, is due to Gutierrez et al. (2002a). P´erez-D´ıaz (2007); Blascoand P´erez-D´ıaz (2019) present a method for computing the singularities of plane curves usinga univariate resultant and characterizing the singularities using its factorization. Notably Rubioet al. (2009) work on rational parametric curves in a ffi ne n -space; they use generalized resultantsto find the parameters of the singular points. Moreover, they characterize the singularities and2ompute their multiplicities.Cox et al. (2011) use the syzygies of the ideal generated by the polynomials that give theparameterization to compute the singularities and their structure. There are state-of-the-art ap-proaches that exploit this idea and relate the problem of computing the singularities with thenotion of the µ -basis of the parametrization, e.g., Jia et al. (2018) and references therein. Chionhand Sederberg (2001) reveal the connection between the implicitization B´ezout matrix and thesingularities of a parametric curve. Bus´e and D’Andrea (2012) present a complete factorizationof the invariant factors of resultant matrices built from birational parameterizations of rationalplane curves in terms of the singular points of the curve and their multiplicity graph. Let usalso mention the important work on matrix methods (Bus´e, 2014; Bus´e and Luu Ba, 2010) forrepresenting the implicit form of parametric curves, that is suitable for numerical computations.Bernardi et al. (2016) use the projection from the rational normal curve to the curve and itsrelation with the secant varieties to the normal curve. Overview of our approach and our contributions.
We introduce
PTOPO , a complete, exact, ande ffi cient algorithm (Alg. 3) for computing the geometric properties and the topology of para-metric curves in R n . Unlike other algorithms, e.g. Alc´azar and D´ıaz-Toca (2010), it makes noassumptions on the input curves, such as the absence of axis-parallel asymptotes, and is applica-ble to any dimension. Nevertheless, it does not handle knots for space curves.If the (proper) parametrization of the curve consists of polynomials of degree d and bitsize τ , then PTOPO outputs a graph isotopic (Boissonnat and Teillaud, 2006, p.184), Alc´azar et al.(2020) to the curve in the embedding space, by performing (cid:101) O B ( nd + nd τ + d ( n + n τ ) + d ( n τ + n ) + n d τ )bit operations in the worst case (Thm. 24). We also provide a Las Vegas variant with expectedcomplexity (cid:101) O B ( d + d ( n + τ ) + d ( n + n τ ) + d ( n τ + n ) + n d τ ) . If n = O (1), the bounds become (cid:101) O B ( N ), where N = max { d , τ } . The vertices of the output graphcorrespond to special points on the curve, in whose neighborhood the topology is not trivial,given by their parameter values. Each edge of the graph is associated with two parameter valuesand corresponds to a unique smooth parametric arc. For an embedding isotopic to the curve, wemap every edge of the abstract graph to the corresponding parametric arc.For plane and space curves, our bound improves the previously known one due to Alc´azar andD´ıaz-Toca (2010) by a factor of (cid:101) O B ( N ). The latter algorithm (Alc´azar and D´ıaz-Toca, 2010)performs some computations in the implicit space. On the contrary, PTOPO is a fundamentallydi ff erent approach since we work exclusively in the parameter space; we do not use a sweep-linealgorithm to construct the isotopic graph. We handle only the parameters that give importantpoints on the curve and thus we avoid performing operations such as univariate root isolation inan extension field or evaluation of a polynomial at an algebraic number.Computing singular points is an essential part of PTOPO (Lem. 18). We chose not to ex-ploit recent methods, e.g., (Blasco and P´erez-D´ıaz, 2019), for this task, to avoid introductingnew variables (T-resultant). We employ older techniques, e.g., (Rubio et al., 2009; P´erez-D´ıaz,2007; Alc´azar and D´ıaz-Toca, 2010), that rely on a bivariate polynomial system, Eq. (2). Wetake advantage of this system’s symmetry and of nearly optimal algorithms for bivariate systemsolving and for computations with real algebraic numbers (Diatta et al., 2018; Bouzidi et al.,2016; Diochnos et al., 2009; Pan and Tsigaridas, 2017). In particular, we introduce an algorithm3or isolating the roots of over-determined bivariate polynomial systems by exploiting the Ratio-nal Univariate Representation (RUR) (Bouzidi et al., 2015b, 2016, 2013) that has worst case andexpected bit complexity that matches the ones for square systems (Thm. 16). These are definitivesteps for obtaining the complexity bounds of Thm. 23 and Thm. 24.Moreover, our bound matches the current state-of-the-art complexity bound, (cid:101) O B ( d + d τ )or (cid:101) O B ( N ), for computing the topology of implicit plane curves Diatta et al. (2018); Kobel andSagralo ff (2014). However, if we want to visualize the graph in 2D or 3D, then we have tocompute a characteristic box (Lem. 20) that contains all the the topological features of the curveand the intersections of the curve with its boundary. In this case, the complexity of PTOPO becomes (cid:101) O B ( N ) (Thm. 23).A preprocessing step of PTOPO consists in finding a proper reparametrization of the curve(if it is not proper). We present explicit bit complexity bounds (Lem. 6) for the algorithm ofP´erez P´erez-D´ıaz (2006) to compute a proper parametrization. Another preprocessing step is toensure that there are no singularities at infinity; Lem. 7 handles this task and provides explicitcomplexity estimates.We additionally consider the case where the embedding of the abstract graph has straight lineedges and not parametric arcs; in particular for plane curves, we show (Cor. 26) that the straightline embedding of the abstract graph to R is already isotopic to the curve. For space curves,the procedure supported by Thm. 28 adds a few extra vertices to the abstract graph, so that thestraight line embedding in R is isotopic to the curve. The extra number of vertices serves inresolving situations where self-crossings occur when continuously deforming the graph to C .We also prove (Thm. 30) that for curves of any dimension we can compute the multiplicities andthe characterization of singular points in the same bit-complexity as computing them (in a LasVegas setting). For that, we use the method by Blasco and P´erez-D´ıaz (2019), which essentiallydoes not require any further computations apart from solving the system that gives the parametersof the singular points (cf.Alc´azar and D´ıaz-Toca (2010)).Last but not least, we provide a certified implementation of PTOPO in maple . The online ver-sion of the implementation handles the topology computation and visualization of plane curves.An updated version handing space curves will be publicly available in the near future.A preliminary version of our work appeared in Katsamaki et al. (2020a). Compared with thisversion, we add all the missing proofs for the algebraic tools that we use in our algorithm, wepresent the isotopic embedding for plane and space curves (Sec. 6), and analyze the complexityof the algorithm of Blasco and P´erez-D´ıaz (2019) to determine the multiplicities and the characterof real singular points (Sec. 7). Organization of the paper.
The next section presents our notation and some useful results neededfor our proofs. In Sect. 3 we give the basic background on rational curves in a ffi ne n -space.We characterize the parametrization by means of injectivity and surjectivity and describe areparametrization algorithm. In Sect. 4 we present the algorithm to compute the singular, ex-treme points, and isolated points on the curve. In Sect. 5 we describe our main algorithm, PTOPO ,that constructs a graph isotopic to the curve in the embedding space and its complexity. In Sect. 6we expatiate on the isotopic embedding for plane and space curves. In Sect. 7, we study the mul-tiplicities and the character of real singular points for curves of arbitrary dimension. Finally, inSect. 8 we give examples and experimental results. https://webusers.imj-prg.fr/~christina.katsamaki/ptopo/ . Notation and Algebraic Tools For a polynomial f ∈ Z [ x ], its infinity norm is equal to the maximum absolute value ofits coe ffi cients. We denote by L ( f ) the logarithm of its infinity norm. We also call the latterthe bitsize of the polynomial. A univariate polynomial is of size ( d , τ ) when its degree is atmost d and has bitsize τ . The bitsize of a rational function is the maximum of the bitsizes ofthe numerator and the denominator. We represent an algebraic number α ∈ C by the isolatinginterval representation . When α ∈ R (resp. C ), it includes a square-free polynomial whichvanishes at α and a (rational) interval (resp. Cartesian products of intervals) containing α and noother root of this polynomial. We denote by O , resp. O B , the arithmetic, resp. bit, complexity andwe use (cid:101) O , resp. (cid:101) O B , to ignore (poly-)logarithmic factors. We denote by res x ( f , g ) the resultantof the polynomials f , g with respect to x . For t ∈ C , we denote by ¯ t its complex conjugate. Weuse [ n ] to signify the set { , . . . , n } .We now present some useful results, needed for our analysis. Lemma 1.
Let A = (cid:80) mi = a i X i , B = (cid:80) ni = b i X i ∈ Z [ X ] of degrees m and n and of bitsizes τ and σ respectively. Let α , . . . , α m be the complex roots of A, counting multiplicities. Then, for any κ = , . . . , m it holds that − m σ − n τ − ( m + n ) log( m + n ) < | B ( α κ ) | < m σ + n τ + ( m + n ) log( m + n ) . Proof.
Following Strzebonski and Tsigaridas (2019), let R = res X ( A ( X ) , Y − B ( X )) ∈ Z [ Y ].Using the Poisson’s formula for the resultant we can write R ( Y ) = a nm (cid:81) m κ = ( Y − B ( α κ )). Themaximum bitsize of the coe ffi cients of R ( Y ) is at most m σ + n τ + ( m + n ) log( m + n ). We observethat the roots of R ( Y ) are B ( α κ ) for κ = , . . . , m . Therefore, using Cauchy’s bound we deducethat 2 − m σ − n τ − ( m + n ) log( m + n ) < | B ( α κ ) | < m σ + n τ + ( m + n ) log( m + n ) . Lemmata 2, 3 restate known results on the gcd computation of various univariate and bivari-ate polynomials.
Lemma 2.
Let f ( X ) , . . . , f n ( X ) ∈ Z [ X ] of sizes ( δ, L ) . We can compute their gcd in worst casecomplexity (cid:101) O B ( n ( δ + δ L )) , or with a Monte Carlo algorithm in (cid:101) O B ( δ + δ L ) , or with a Las Vegasalgorithm in (cid:101) O B ( n ( δ + δ L )) .Proof. These are known results von zur Gathen and Gerhard (2013); we repeat the argumentsadapted to our notation.
Worst case:
We compute g by performing n consecutive gcd computations, that isgcd( f , gcd( f , gcd( · · · , gcd( f n − , f n ))) . Since each gcd computation costs (cid:101) O B ( δ + δ L ) (Bouzidi et al., 2015b, Lem.4), the result for thiscase follows. Monte Carlo:
We perform one gcd computation by allowing randomization. If we chooseintegers a , . . . , a n independently at random from the set { , . . . , Kd } , where K = O (1), we getthat gcd( f , . . . , f n ) = gcd( f , f + a f + · · · + a n f n ) in Z [ x ], with probability ≥ / Q . To computethe gcd in Z we need to multiply with the gcd of the leading coe ffi cients of f , f + a f + · · · + a n f n ffi cient since the leadingcoe ffi cient of the gcd in Z [ x ] divides the leading coe ffi cients of the two polynomials and by(von zur Gathen and Gerhard, 2013, Cor. 6.10) the monic gcd of two polynomials in Q [ x ] is equalto their gcd in Z [ x ] divided by their leading coe ffi cient. The gcd of the two leading coe ffi cientsof f , f + a f + · · · + a n f n is an integer of bitsize (cid:101) O ( L ), therefore this does not pollute the totalcomplexity.We compute g ∗ = gcd( f , f + a f + · · · + a n f n ). Notice that the polynomial f + a f + · · · + a n f n is asymptotically of size ( δ, L ). So, it takes (cid:101) O B ( δ + δ L ) to find g ∗ , using the probabilistic algorithmin Sch¨onhage (1988). Las Vegas:
We can reduce the probability of failure in the Monte Carlo variant of the gcdcomputation to zero, by performing n exact divisions. In particular, we check if g ∗ divides f , . . . , f n . Using (von zur Gathen and Gerhard, 2013, Ex.10.21), the bit complexity of theseoperations is in total (cid:101) O B ( n ( δ + δτ )). Lemma 3.
Let f ( X , Y ) , . . . , f n ( X , Y ) ∈ Z [ X , Y ] of bidegrees ( δ, δ ) and L ( f i ) = L. We can com-pute their gcd in worst case complexity (cid:101) O B ( n ( δ + δ L )) , or with a Monte Carlo algorithm in (cid:101) O B ( δ + δ L ) , or with a Las Vegas algorithm in (cid:101) O B ( n ( δ + δ L )) .Proof. The straightforward approach is to perform n consecutive gcd computations, that isgcd( f , gcd( f , gcd( · · · , gcd( f n − , f n ))) . To accelerate the practical complexity we should sort f i in increasing order with respect to theirdegree. Each gcd computation costs (cid:101) O B ( δ + δ L ) (Bouzidi et al., 2016, Lem. 5), so the totalworst case cost is (cid:101) O B ( n δ + n δ L ).Alternatively, we consider the operation gcd( f , (cid:80) nk = a k f k ), where a k are random integers,following (von zur Gathen and Gerhard, 2013, Thm. 6.46). The expected cost of this gcd is (cid:101) O B ( δ + δ L ). To see this, notice that we can perform a bivariate gcd in expected time (cid:101) O ( δ )(von zur Gathen and Gerhard, 2013, Cor. 11.12), over a finite field with enough elements, andthe bitsize of the result is (cid:101) O ( δ + L ) Mahler (1962).Then, for a Las Vegas algorithm, using exact division, we test if the resulting polynomialdivides all f i , for 2 ≤ i ≤ n . This costs (cid:101) O B ( n ( δ + δ L )), by adapting (von zur Gathen andGerhard, 2013, Ex.10.21) to the bivariate case.
3. Rational curves
Following closely Alc´azar and D´ıaz-Toca (2010), we introduce basic notions for rationalcurves. Let (cid:101) C be an algebraic curve over C n , parametrized by the map φ : C (cid:57)(cid:57)(cid:75) (cid:101) C t (cid:55)→ (cid:0) φ ( t ) , . . . , φ n ( t ) (cid:1) = (cid:16) p ( t ) q ( t ) , . . . , p n ( t ) q n ( t ) (cid:17) , (1)where p i , q i ∈ Z [ t ] are of size ( d , τ ) for i ∈ [ n ], and (cid:101) C is the Zariski closure of Im( φ ). We call φ ( t )a parametrization of (cid:101) C .We study the real trace of ˜ C , that is C : = ˜ C ∩ R n . A parametrization φ is chatacterized bymeans of properness (Sec. 3.1) and normality (Sec. 3.2). To ensure these properties, one can6eparametrize the curve, i.e., apply a rational change of parameter to the given parametrization.We refer to (Sendra et al., 2008, Ch. 6) for more details on reparameterization.Without loss of generality, we assume that no component of the parametrization φ is constant;otherwise we could embed ˜ C in a lower dimensional space. We consider that φ is in reduced form ,i.e., gcd( p i ( t ) , q i ( t )) =
1, for all i ∈ [ n ]. The point at infinity, p ∞ , is the point on C we obtainfor t → ±∞ (if it exists). For a parametrization φ , we consider the following system of bivariatepolynomials: h i ( s , t ) = p i ( s ) q i ( t ) − q i ( s ) p i ( t ) s − t , for i ∈ [ n ] . (2) Remark 4.
The h i ’s are polynomials since ( s , s ) is a root of the numerator for every s. Also,h i ( t , t ) = φ (cid:48) i ( t ) q i ( t ) for i ∈ [ n ] (Gutierrez et al., 2002b, Lem. 1.7).3.1. Proper parametrization A parametrization is proper if φ ( t ) is injective for almost all points on (cid:101) C . In other words,almost every point on (cid:101) C is the image of exactly one parameter value (real or complex). Forother equivalent definitions of properness we refer to (Sendra et al., 2008, Ch. 4), Rubio et al.(2009). The following condition (Alc´azar and D´ıaz-Toca, 2010, Thm. 1) leads to an algorithm forchecking properness: a parametrization is proper if and only if deg(gcd( h ( s , t ) , . . . , h n ( s , t ))) = Lemma 5.
There is an algorithm that checks if a parametrization φ is proper in worst-case bitcomplexity (cid:101) O B ( n ( d + d τ )) and in expected bit complexity (cid:101) O B ( n ( d + d τ )) .Proof. The construction of all h i costs O B ( nd τ ). We need to check if deg(gcd( h ( s , t ) , . . . , h n ( s , t ))) = φ is a not a proper parametrization, then there always exists a parametrization ψ ∈ Z ( t ) n and R ( t ) ∈ Z ( t ) such that ψ ( R ( t )) = φ ( t ) and ψ is proper (Sendra et al., 2008, Thm. 7.6). Thereare various algorithms for obtaining a proper parametrization, e.g., Sederberg (1986); Gutierrezet al. (2002a); Sendra et al. (2008); P´erez-D´ıaz (2006); Gao and Chou (1992). We consider thealgorithm in P´erez-D´ıaz (2006) for its simplicity; its pseudo-code is in Alg. 1. Lemma 6.
Consider a non-proper parametrization of a curve C , consisting of univariate poly-nomials of size ( d , τ ) . Alg. 1 computes a proper parametrization of C , involving polynomials ofdegree at most d and bitsize O ( d + d τ ) , in (cid:101) O B ( n ( d + d τ )) , in the worst case.Proof. The algorithm first computes the bivariate polynomials H , . . . H n . They have bi-degreeat most ( d , d ) and bitsize at most 2 τ +
1. Then, we compute their gcd, that we denote by H , in (cid:101) O B ( n ( d + d τ )) (Lem. 3). By Mahler (1962) and (Basu et al., 2003, Prop. 10.12) we have that L ( H ) = O ( d + τ ), which is also the case for C j ( s ).If the degree of H is one, then the parametrization is already proper and we have nothing todo. Otherwise, we consider H as a univariate polynomial in s and we find two of its coe ffi cientsthat are relatively prime, using exact division. The complexity of this operation is m · (cid:101) O B ( d + d τ ) = (cid:101) O B ( d + d τ ) (von zur Gathen and Gerhard, 2013, Ex. 10.21).Subsequently, we perform n resultant computations to get L , . . . L n . From these we obtainthe rational functions of the new parametrization. We focus on the computation of L . The same7rguments hold for all L i . The bi-degree of L ( s , x ) is ( d , d ) (Basu et al., 2003, Prop. 8.49) and L ( L ) = O ( d + d τ ) (Basu et al., 2003, Prop. 8.50); the latter dictates the bitsize of the newparametrization.To compute L , we consider F and G as univariate polynomials in t and we apply a fastalgorithm for computing the univariate resultant based on subresultants Lickteig and Roy (2001);it performs (cid:101) O ( d ) operations. Each operation consists of multiplying bivariate polynomials of bi-degree ( d , d ) and bitsize O ( d + d τ ); so it costs (cid:101) O B ( d + d τ ). We compute the resultant in (cid:101) O B ( d + d τ ). We multiply the latter bound by n to conclude the proof. Algorithm 1:
Make Proper ( φ ) Input:
A parametrization φ ∈ Z ( t ) n as in Eq. (1) Output:
A proper parametrization ψ = ( ψ , . . . , ψ n ) ∈ Z ( t ) n for i ∈ [ n ] do H i ( s , t ) ← p i ( s ) q i ( t ) − p i ( t ) q i ( s ) ∈ Z [ s , t ] ; H ← gcd( H , . . . , H n ) = C m ( t ) s m + · · · + C ( t ) ∈ ( Z [ t ])[ s ] if m = then return φ ( t ) ; Find k , l ∈ [ m ] such that:deg(gcd( C k ( t ) , C l ( t ))) = C k ( t ) C l ( t ) (cid:60) Q R ( t ) ← C k ( t ) C l ( t ) r ← deg( R ) = max { deg( C k ) , deg( C l ) } G ← s C l ( t ) − C k ( t ) for i ∈ [ n ] do F i ← xq i ( t ) − p i ( t ) L i ( s , x ) ← res t ( F i ( t , x ) , G ( t , s )) = ( ˜ q i ( s ) x − ˜ p i ( s )) r return ψ ( t ) = (cid:0) ˜ p ( t )˜ q ( t ) , . . . , ˜ p n ( t )˜ q n ( t ) (cid:1) Normality of the parametrization concerns the surjectivity of the map φ . The parametrization φ ( t ) is R -normal if for all points p on C there exists t ∈ R such that φ ( t ) = p . When theparametrization is not R -normal, the points that are not in the image of φ for t ∈ R are p ∞ (if itexists) and the isolated points that we obtain for complex values of t (Recio, 2007, Prop. 4.2).An R -normal reparametrization does not always exist. We refer to (Sendra et al., 2008, Sect. 7.3)for further details.However, if p ∞ exists, then we reparametrize the curve to avoid possible singularities atinfinity. The point p ∞ exists if deg( p i ) ≤ deg( q i ), for all i ∈ [ n ]. Lemma 7. If p ∞ exists, then we can reparametrize the curve with a linear function to ensurethat p ∞ is not a singular point, using a Las Vegas algorithm in expected time (cid:101) O B ( n ( d + d τ )) . Thenew parametrization involves polynomials of size ( d , (cid:101) O ( d + τ )) .Proof. The point at infinity depends on the parametrization. So, for this proof, let us denote thepoint at infinity of φ by p φ ∞ .The reparametrization consists in choosing t ∈ R and applying the map r : t (cid:55)→ t t + t − t to φ ,to obtain a new parametrization, ψ = φ ◦ r . The point at infinity of the new parametrization is8 ψ ∞ = φ ( t ). We need to ensure that p ψ ∞ = φ ( t ) is not singular. There are O ( d ) singular points,so we choose t uniformly at random from the set { , . . . , Kd } where K = O (1). Then, withprobability ≥ / φ ( t ) is not singular and p ψ ∞ is also not singular. The bound on the possiblevalues of t implies that the bitsize of t is O (lg( d )).We compute the new parametrization, ψ , in (cid:101) O B ( n ( d + d τ )) using multipoint evaluation andinterpolation, by exploiting the fact that the polynomials in ψ have degrees at most d and bitsize (cid:101) O ( d + τ ).For a Las Vegas algorithm we need to check if φ ( t ) is a cusp or a multiple point. For theformer, we evaluate φ (cid:48) at t (see Rem. 4). This costs (cid:101) O B ( nd τ ) (Bouzidi et al., 2013, Lem. 3). Forthe latter, we check if deg(gcd( φ ( t ) q ( t ) − p ( t ) , . . . , φ ( t ) q ( t ) − p ( t ))) = (cid:101) O B ( n ( d + d τ ))(Lem. 2). If φ (cid:48) ( t ) is not the zero vector and the degree of the gcd is zero, then φ ( t ) is notsingular. Remark 8.
Since the reparametrizing function in the previous lemma is linear, it does not a ff ectproperness (Sendra et al., 2008, Thm. 6.3).
4. Special points on the curve
We consider a parametrization φ of C as in Eq. (1), such that φ is proper and there are nosingularities at infinity. We highlight the necessity of these assumptions when needed. We detectthe parameters that generate the special points of C , namely the singular, the isolated, and theextreme points. We identify the values of the parameter for which φ is not defined, namely thepoles (see Def. 9); in presence of poles, C consists of multiple components. Definition 9.
The parameters for which φ ( t ) is not defined are the poles of φ . The sets of polesover the complex and the reals are: T C P = { t ∈ C : (cid:89) i ∈ [ n ] q i ( t ) = } and T R P = T C P ∩ R . We consider the solution set S of system (2) over C : S = { ( t , s ) ∈ C : h i ( t , s ) = i ∈ [ n ] } . Remark 10.
Notice that when φ is in reduced form, if ( s , t ) ∈ S and ( s , t ) ∈ ( C \ T C P ) × C , thenalso t (cid:60) T C P (Rubio et al., 2009, (in the proof of) Lem. 9). Next, we present some well known results Rubio et al. (2009); Sendra et al. (2008) that weadapt in our notation.
Singular points . Quoting Manocha and Canny (1992), ”Algebraically, singular points are pointson the curve, in whose neighborhood the curve cannot be represented as an one-to-one and C ∞ bijective map with an open interval on the real line”. Geometrically, singularities correspondto shape features that are known as cusps and self-intersections of smooth branches. Cusps arepoints on the curve where the tangent vector is the zero vector. This is a necessary and su ffi cientcondition when the parametrization is proper Manocha and Canny (1992). Self-intersections are multiple points, i.e., points on C with more than one preimages.9 emma 11. The set of parameters corresponding to real cusps is T C = (cid:110) t ∈ R \ T R P : ( t , t ) ∈ S (cid:111) . The set of parameters corresponding to real multiple points is T M = { t ∈ R \ T R P : ∃ s (cid:44) t , s ∈ R such that ( t , s ) ∈ S } . Proof.
The description of T C is an immediate consequence of Rem. 4. It states that h i ( t , t ) = φ (cid:48) i ( t ) q i ( t ), for i ∈ [ n ].Now let p = φ ( t ) be a multiple point on C . Then, there is s ∈ R \ T R P with φ ( t ) = φ ( s ) ⇒ h i ( t , s ) = i ∈ [ n ] and so t ∈ T M . Conversely, let t ∈ T M and s (cid:44) t , s ∈ R suchthat h i ( t , s ) = i ∈ [ n ]. From (Rubio et al., 2009, (in the proof of) Lem. 9), when φ is inreduced form, if ( t , s ) ∈ S and ( t , s ) ∈ ( R \ T R P ) × R , then also s (cid:60) T R P . So, h i ( t , s ) = ⇔ p i ( t ) q i ( t ) = p i ( s ) q i ( s ) for all i ∈ [ n ], and thus p = φ ( t ) = φ ( s ) is real multiple point.Notice that T C and T M are not necessarily disjoint, for at the same point we may have both cuspsand smooth branches that intersect. Isolated points . An isolated point on a real curve can only occur for complex values of theparameter. The point at infinity is not isolated because it is the limit of a sequence of real points.
Lemma 12.
The set of parameters generating isolated points of C is T I = { t ∈ C \ ( R ∪ T C P ) : ( t , t ) ∈ S and (cid:64) s ∈ R s.t. ( t , s ) ∈ S } . Proof.
Let p = φ ( t ) ∈ R n be an isolated point, where t ∈ C \ ( R ∪ T C P ). Notice that p is also amultiple point, since it holds that φ i ( t ) = φ i ( t ) = φ i ( t ) for i ∈ [ n ]. Thus, h i ( t , t ) = i ∈ [ n ]and ( t , t ) ∈ S . Moreover, since p is isolated, there are no real branches through p and there doesnot exist s ∈ R such that φ ( t ) = φ ( s ) ⇒ h i ( t , s ) =
0, for all i ∈ [ n ]. So, t ∈ T I .Conversely, let ( t , t ) ∈ S with t ∈ C \ R ∪ T C P . Since φ is in reduced form, we have that t (cid:60) P C (Rubio et al., 2009, (in the proof of) Lem. 9), therefore h i ( t , t ) =
0, for all i ∈ [ n ], impliesthat φ ( t ) = φ ( t ) = φ ( t ) ∈ R n . Since there does not exist s ∈ R with φ ( t ) = φ ( s ), p is an isolatedpoint on C . Extreme points . Consider a vector (cid:126)δ and a point on C whose tangent vector is parallel to (cid:126)δ . If thepoint is not singular, then it is an extreme point of C with respect to (cid:126)δ . We compute the extremepoints with respect to the direction of each coordinate axis. Rem. 4 leads to the following lemma: Lemma 13.
The set of parameters generating extreme points is T E = (cid:8) t ∈ R \ T R P : (cid:89) i ∈ [ n ] h i ( t , t ) = and t (cid:60) T C ∪ T M (cid:9) . lgorithm 2: Special Points ( φ ) Input:
Proper parametrization φ ∈ Z ( t ) n without singularity at infinity, as in Eq. (1) Output:
Real poles and parameters that give real cusps, multiple, isolated and extremepoints. /* The subroutines SOLVE R and SOLVE C return the solution set of aunivariate polynomial or a system of polynomials over the realand complex numbers resp. */ Compute polynomials h ( s , t ) , . . . , h n ( s , t ) T R P ← (cid:83) i ∈ n SOLVE R ( q i ( t ) = T C P ← (cid:83) i ∈ [ n ] SOLVE C ( q i ( t ) = S ← SOLVE C ( h ( s , t ) = , . . . , h n ( s , t ) = T C , T M , T I , W ← ∅ for ( s , t ) ∈ S do if s = t and s ∈ R \ T R P then T C ← T C ∪ { t } else if s (cid:44) t then if s ∈ R \ T R P then if t ∈ R then T M ← T M ∪ { t } else W ← W ∪ { t } else if s = t and s (cid:60) T C P then T I ← T I ∪ { t } T I ← T I \ W /* Extreme points */ T E ← (cid:83) i ∈ n SOLVE R ( h i ( t , t ) = T E ← T E \ ( T E ∩ ( T C ∪ T M )) From Lemmata 11, 12, and 13, it follows that given a proper parametrization φ withoutsingular points at infinity, we can easily find the poles and the set of parameters generatingcusps, multiple, extreme, and isolated points. We do so, by solving an over-determined bivariatepolynomial system and univariate polynomial equations. Then, we classify the parameters thatappear in the solutions, by exploiting the fact the system is symmetric. For sake of completeness,we describe the procedure in Alg. 2.To compute the RUR of an overdetermined bivariate system (Thm. 16), we employ Lem. 14and Prop. 15, which adapt the techniques used in Bouzidi et al. (2016) to our setting. Lemma 14.
Let f , g , h , . . . , h n ∈ Z [ X , Y ] with degrees bounded by δ and bitsize of coe ffi cientsbounded by L. Computing a common separating element in the form X + α Y , α ∈ Z for then + { f = g = } , { f = h i = } , i = . . . n needs (cid:101) O B ( n ( δ + δ L )) bit operations in the worst case, and (cid:101) O B ( n ( δ + δ L )) in the expected case witha Las Vegas Algorithm. Moreover, the bitsize of α does not exceed log(2 n δ ) . roof. A straightforward strategy consists in running simultaneously Algorithm 5 (worst case) orAlgorithm 5’ (Las Vegas) from Bouzidi et al. (2016) on all the systems. The only modificationsneeded are that the values of α to be considered are less than 2 n δ (twice a bound on the totalnumber of solutions of all the systems) and that the exit test is valid if and only if it is valid forall the systems. Proposition 15.
Let f , g ∈ Z [ X , Y ] with degrees bounded by δ and coe ffi cients’ bitsizes boundedby L. We can compute a rational parameterization { h ( T ) , X = h X ( T ) h ( T ) , Y = h Y ( T ) h ( T ) } of f , g withh , h , h X , h Y ∈ Z [ T ] with degrees less than δ and coe ffi cients’ bitsizes in (cid:101) O ( δ ( L + δ )) , in (cid:101) O B ( δ ( L + δ )) bit operations in the worst case and (cid:101) O B ( δ ( L + δ )) expected bit operations with a Las VegasAlgorithm.Proof. Algorithms 6 and 6’ from Bouzidi et al. (2016) compute a RUR decomposition of f = g = (cid:101) O B ( δ ( L + δ )) bit operations in the worst case and (cid:101) O B ( δ ( L + δ )) expected bit operationswith a Las Vegas Algorithm respectively. They provide s (cid:54) δ parameterizations in the form { h i ( T ) , h i , X ( T ) h i , ( T ) , h i , Y ( T ) h i , ( T ) } , where i = .. s , with the following properties: • (cid:81) si = h i is a polynomial of degree at most δ with coe ffi cients of bitsize (cid:101) O ( δ L + δ ). • The degrees of h i , ( T ) , h X , ( T ) and h Y , ( T ) are less then the degree of h i . • The coe ffi cients’ bitsizes of h i , ( T ) , h X , ( T ) and h Y , ( T ) are in (cid:101) O B ( δ L + δ ).Also, s (cid:89) i = h i , (cid:80) nn = h j , X (cid:81) i (cid:44) j h i (cid:80) nn = h j , (cid:81) i (cid:44) j h i , (cid:80) nn = h j , Y (cid:81) i (cid:44) j h i (cid:80) nn = h j , (cid:81) i (cid:44) j h i is a rational parameterization of the system { f = g = } , defined by polynomials of degree lessthan δ with coe ffi cients of bitsizes (cid:101) O ( δ ( L + δ )) and can be computed from the RUR decompositionperforming O ( s ) multiplications of polynomials of degree at most δ with coe ffi cients of bitsize (cid:101) O ( δ ( L + δ )), which requires (cid:101) O B ( δ ( L + δ )) bit operations. Theorem 16.
There exists an algorithm that computes the RUR and the isolating boxes of theroots of the system { h ( s , t ) = · · · = h n ( s , t ) = } with worst-case bit complexity (cid:101) O B ( n ( d + d τ )) .There is also a Las Vegas variant with expected complexity (cid:101) O B ( d + nd + d τ + nd τ ) .Proof. Assume that we know a common separating linear element (cid:96) ( s , t ) = (cid:96) + (cid:96) s + (cid:96) t thatseparates the roots of the n -1 systems of bivariate polynomial equations { h = h = } , { h = h i = } , for 3 ≤ i ≤ n . We can compute (cid:96) with (cid:101) O B ( n ( d + d τ )) bit operations in the worst caseand with (cid:101) O B ( n ( d + d τ )) expected bit operations with a Las Vegas algorithm (Lem. 14).We denote by { r ( T ) , r s ( T ) r I ( T ) , r t ( T ) r I ( T ) } a RUR for { h = h = } with respect to (cid:96) . In addition, for i = . . . n , let { r i ( T ) , r i , s ( T ) r i , I ( T ) , r i , t ( T ) r i , I ( T ) } be the RUR of { h = h i = } , also with respect to (cid:96) . We cancompute all these representations with (cid:101) O B ( n ( d + d τ )) bit operations in the worst case, and with (cid:101) O B ( n ( d + d τ )) in expected case with a Las Vegas algorithm (Lem. 15).Then, for the system { h = h = . . . = h n = } we can define a rational parameterization { χ ( T ) , r s ( T ) r I ( T ) , r t ( T ) r I ( T ) } , where χ ( T ) = gcd( r ( T ) , r ( T ) , . . . , r n ( T ) , r s ( T ) r , I ( T ) − r , s ( T ) r I ( T ) , r t ( T ) r , I ( T ) − r , t ( T ) r I ( T ) ,... r s ( T ) r n , I ( T ) − r n , s ( T ) r I ( T ) , r t ( T ) r n , I ( T ) − r n , t ( T ) r I ( T )) . n − d and coe ffi cients of bitsizes in (cid:101) O ( d τ ) which needs (cid:101) O B ( n ( d + d τ )) bit operations in the worst case. Isolating the roots of such a parameterization requires (cid:101) O B ( d + d τ ) according to Alg. 7 from Bouzidi et al. (2016). Remark 17 (RUR and isolating interval representation) . If we use Thm.16 to solve the over-determined bivariate system of the h i polynomials of Eq. (2), then we obtain in the output a RURfor the roots, which is as follows: There is a polynomial χ ( T ) ∈ Z [ T ] of size ( O ( d ) , (cid:101) O ( d + d τ )) and a mapping V ( χ ) → V ( h , . . . , h n ) T (cid:55)→ (cid:0) r s ( T ) r I ( T ) , r t ( T ) r I ( T ) (cid:1) , (3) that defines an one-to-one correspondence between the roots of χ and those of the system. Thepolynomials r s , r t , and r I are in Z [ T ] and have also size ( O ( d ) , (cid:101) O ( d + d τ )) .Taking into account the cost to compute this parametrization of the solutions (Thm.16), wecan also compute at no extra cost the resultant of { h , h } with respect to s or t. Notice that bothresultants are the same polynomial, since the system is symmetric. Let R s ( t ) = res s ( h , h ) . It isof size ( O ( d ) , O ( d + d τ )) (Basu et al., 2003, Prop. 8.46).Under the same bit complexity, we can su ffi ciently refine the isolating boxes of the solutionsof the bivariate system (computed in Thm.16), so that every root ( r s ( ξ ) r I ( ξ ) , r t ( ξ ) r I ( ξ ) ) , where χ ( ξ ) = , hasa representation as a pair of algebraic numbers in isolating interval representation: (( R s , I ,ξ × I ,ξ ) , ( R s , J ,ξ × J ,ξ )) . (4) Both coordinates are roots of the same polynomial. Moreover, I ,ξ , J ,ξ are empty sets when thecorresponding algebraic number is real. Therefore, we can immediately distinguish between realand complex parameters. At the same time, we associate to each isolating box of a root of R s thealgebraic numbers ρ = ( χ, I ρ × J ρ ) for whom it holds that r s ( ρ ) r I ( ρ ) projects inside this isolating box.We can interchange between the two of representations in constant time, and this will simplifyour computations in the sequel. Lemma 18.
Let C be a curve with a proper parametrization φ ( t ) as in Eq. (1), that has nosingularities at infinity. We compute the real poles of φ and the parameters corresponding tosingular, extreme, and isolated points of C in worst-case bit complexity (cid:101) O B ( nd + nd τ + d ( n + n τ ) + d ( n τ + n ) + n d τ ) , and using a Las Vegas algorithm in expected bit complexity (cid:101) O B ( d + d ( n + τ ) + d ( n + n τ ) + d ( n τ + n ) + n d τ ) . Proof.
The proof is an immediate consequence of the following: • We compute all h i ∈ Z [ s , t ] in (cid:101) O B ( nd τ ): To construct each h i we perform d multiplicationsof numbers of bitsize τ ; the cost for this is (cid:101) O B ( d τ ). The bi-degree of each is at most ( d , d ) and L ( h i ) ≤ τ + = O ( τ ). 13 The real poles of φ are computed in (cid:101) O B ( n ( d + d τ )) : To find the poles of φ , we isolate thereal roots of each polynomial q i ( t ), for i ∈ [ n ]. This costs (cid:101) O B ( n ( d + d τ )) Pan and Tsigaridas(2017). Then we sort the roots in (cid:101) O B ( n d n ( d + d τ )) = (cid:101) O B ( n ( d + d τ )). • The parameters corresponding to cusps, multiple and isolated points of C are computed in (cid:101) O B ( n ( d + d τ )) : We solve the bivariate system (2) in (cid:101) O B ( n ( d + d τ )) or in expected time (cid:101) O B ( d + nd + d τ + nd τ )(Thm. 16). Then we have a parametrization of the solutions of the bivariate system (2) of theform (3) and in the same time of the form (4) (see Rem. 17). Some solutions ( s , t ) ∈ S may notcorrespond to points on the curve, since s , t can be poles of φ . Notice that from Rem. 10, s and t are either both poles or none of them is a pole. We compute g s = gcd( R s , Q ), where Q ( t ) = (cid:81) i ∈ [ n ] q i ( t ), and the gcd-free part of R s with respect to Q . This is done in (cid:101) O B (max { n , d } ( nd τ + nd τ )) (Bouzidi et al., 2015b, Lem. 5).Every root of R ∗ s is an algebraic number of the form ( R s , I ,ξ × I ,ξ ), for some ξ that is root of χ .We can easily determine if it corresponds to a cusp, a multiple or an isolated point; when real(i.e., I ,ξ = ∅ ) it corresponds to a cusp of C if and only if (( R s , I ,ξ ) , ( R s , I ,ξ )) is in S . Otherwise,it corresponds to a multiple point. When it is complex (i.e., I ,ξ (cid:44) ∅ ), it corresponds to an isolatedpont of C if and only if (( R s , I ,ξ × I ,ξ ) , ( R s , I ,ξ × ( − I ,ξ ))) ∈ S and there is no root in S of theform (( R s , I ,ξ × I ,ξ ) , ( R s , J ,ξ (cid:48) )). • The parameters corresponding to extreme points of C are computed in (cid:101) O B ( d n τ + d ( n τ + n ) + d n τ ) : For all i ∈ [ n ], h i ( t , t ) is a univariate polynomial of size ( O ( d ) , O ( τ )). Then, H ( t ) = (cid:81) i ∈ [ n ] h i ( t , t ) isof size ( O ( nd ) , (cid:101) O ( n τ )). The parameters that correspond to extreme points are among the roots of H ( t ). To make sure that poles and parameters that give singular points are excluded, we computegcd( H , Q · R s ), where Q ( t ) = (cid:81) i ∈ [ n ] q i ( t ), and the gcd-free part of H with respect to Q · R s , say H ∗ .Since Q · R s is a polynomial of size ( d + nd , ( d + n ) τ ), the computation of the gcd and the gcd-freepart costs (cid:101) O B ( n ( d τ + nd τ + n d τ )) (Bouzidi et al., 2015b, Lem. 5). Then, H = gcd( H , Q · R s ) H ∗ ,and the real roots of H ∗ give the parameters that correspond to extreme points. We isolate thereal roots of H ∗ in (cid:101) O B ( n ( d + d τ )), since it is a polynomial of size ( O ( nd ) , (cid:101) O ( n ( d + τ ))). PTOPO : Topology and Complexity
We present
PTOPO , an algorithm to construct an abstract graph G that is isotopic (Boissonnatand Teillaud, 2006, p.184) to C when we embed it in R n . We emphasize that, currently, we donot treat / compute knots in the case of space curves. The embedding consists of a graph whosevertices are points on the curve given by their parameter values. The edges are smooth parametricarcs that we can continuously deform to branches of C without any topological changes. Weneed to specify a bounding box in R n inside which the constructed graph results in an isotopicembedding to C . We comment at the end of the section on the case where an arbitrary box isprovided at the input. We determine a bounding box in R n , which we call characteristic , thatcaptures all the topological information of C : Definition 19.
A characteristic box of C is a box enclosing a subset of R n that intersects allcomponents of C and contains all its singular, extreme, and isolated points. Let B C be a characteristic box of C . If C is bounded, then C ⊂ B C . If C is unbounded, thenthe branches of C that extend to infinity intersect the boundary of B C . A branch of the curveextends to infinity if for t → t , it holds || φ ( t ) || > M , for any M >
0, where t ∈ R ∪ {∞} .14em. 20 computes a characteristic box using the degree and bitsize of the polynomials in theparametrization of Eq. (1). Lemma 20.
Let C be a curve with a parametrization as in Eq. (1). For b = d ( τ + log d ) = O ( d τ ) , B C = [ − b , b ] n is a characteristic box of C .Proof. We estimate the maximum and minimum values of φ i , i ∈ [ n ], when we evaluate it at theparameter values that correspond to special points and also at each pole that is not a root of q i .Let t be a parameter that corresponds to a cusp or an extreme point with respect to the i -th direction. Then, it is a root of φ (cid:48) i ( t ). Let N ( t ) = p (cid:48) i ( t ) q i ( t ) − p i ( t ) q (cid:48) i ( t ) the numerator of φ (cid:48) i ( t ). Then N ( t ) =
0. The degree of N ( t ) is ≤ d − L ( N ) ≤ τ + log d + . From Lem. 1we conclude that | p i ( t ) | ≤ d τ + d log( d ) + (3 d −
1) log(3 d − + d − τ . Analogously, it holds that | q i ( t ) | ≥ − d τ − d log( d ) − (3 d −
1) log(3 d − − d + τ . Therefore, | φ i ( t ) | ≤ d τ + d log( d ) + (3 d −
1) log(3 d − + d − τ ) . Now, let ( t , t ) be two parameters corresponding to a multiple point of C , i.e., ( t , t ) is a root ofthe bivariate system in Eq. (2). Take any j , k ∈ [ n ] with j (cid:44) k and let R ( t ) = res s ( h j , h k ). It holdsthat R ( t ) =
0. The degree of R is ≤ d and L ( R ) ≤ d ( τ + log( d ) + log( d + +
1) (Basu et al.,2003, Prop. 8.29). Applying Lem. 1, we deduce that | φ i ( t ) | ≤ d ( τ + log( d ) + log( d + + + d τ + (2 d + d ) log(2 d + d ) . Let t be a pole of φ with q j ( t ) =
0, for some j (cid:44) i . If φ i ( t ) is defined, applying Lem. 1 gives | φ i ( t ) | ≤ d τ + d log 2 d . To conclude, we take the maximum of the three bounds. However, to simplify notation, weslightly overestimate the latter bound.The vertices of the embedded graph must include the singular and the isolated points of C .Additionally, to rigorously visualize the geometry of C , we consider as vertices the extremepoints of C , with respect to all coordinate directions, as well as the intersections of C with theboundary of the bounding box. We label the vertices of G using the corresponding parametervalues generating these points and we connect them accordingly. Alg. 3 presents the pseudo-code of PTOPO and here we give some more details on the various steps. We construct G asfollows:First, we compute the poles and the sets T C , T M , T E , and T I of parameters corresponding to“special points”. Then, we compute the characteristic box of C , say B C . We compute the set T B of parameters corresponding to the intersections of C with the boundary of B C (if any). Lem. 21describes this procedure and its complexity. Lemma 21.
Let B = [ l , r ] ×· · ·× [ l n , r n ] in R n and L ( l i ) = L ( r i ) = σ , for i ∈ [ n ] . We can find theparameters that give the intersection points of φ with the boundary of B in (cid:101) O B ( n d + n d ( τ + σ )) .Proof. For each i ∈ [ n ] the polynomials q i ( t ) l i − p i ( t ) = q i ( t ) r i − p i ( t ) = O ( d ) , O ( τ + σ )). So, we compute isolating intervals for all their real solutions in (cid:101) O B ( d ( τ + σ )) Pan(2002). For any root t of each of these polynomials, since φ is in reduced form (by assumption),we have t (cid:60) T R P . We check if φ j ( t ) ∈ [ l j , r j ], j ∈ [ n ] \ i . This requires 3 sign evaluations ofunivariate polynomials of size ( d , τ + σ ) at all roots of a polynomial of size ( d , τ + σ ). The bit15omplexity of performing these operations for all the roots is (cid:101) O B ( d + d ( τ + σ )) (Strzebonskiand Tsigaridas, 2019, Prop. 6). Since we repeat this procedure n − i ∈ [ n ], thetotal cost is (cid:101) O B ( n d + n d ( τ + σ )).We partition T C ∪ T M ∪ T E ∪ T I ∪ T B into groups of parameters that correspond to the samepoint on C . For each group, we add a vertex to G if and only if the corresponding point is insidethe bounding box B ; for the characteristic box it is inside by construction. Lemma 22.
The graph G has κ = O ( d + nd ) vertices, which can be computed using O ( d + nd ) arithmetic operations.Proof. Since T B ∩ T M = ∅ and T E ∩ T M = ∅ , to each parameter in T B and T E corresponds a uniquepoint on C . So for every t ∈ T B ∪ T E we add a vertex to G , labeled by the respective parameter.Next, we group the parameters in T C ∪ T M ∪ T I that give the same point on C and we add for eachgroup a vertex at G labeled by the corresponding parameter values.Grouping of the parameters is done as follows: For every t ∈ T C ∪ T M we add a vertex to G labeled by the set { s ∈ R : ( s , t ) ∈ S } ∪ { t } and for every t ∈ T I we add a vertex to G labeled bythe set { s ∈ C : ( s , t ) ∈ S } ∪ { t } . Notice that we took into account Rem.10. We compute these setssimply by reading the elements of S .It holds that T B = O ( nd ), T E = O ( nd ) and | S | = O ( d ). Since for each vertex, we can find theparameters that give the same point in constant time, the result follows.We denote by v , . . . , v κ the vertices (with distinct labels) of G and by λ ( v ) , . . . , λ ( v κ ) theirlabel sets (i.e., the parameters that correspond to each vertex). Let T be the sorted list of param-eters in T C ∪ T M ∪ T E ∪ T B . If for two consecutive elements t < t in T , there exists a pole s ∈ T R P such that t < s < t , then we split T into two lists: T containing the elements ≤ t and T containing the elements ≥ t . We continue recursively for T and T , until there are no polesbetween any two elements of the resulting list. This procedure partitions T into T , . . . , T (cid:96) .To add edges to G , we consider each T i with more than one element, where i ∈ [ (cid:96) ], indepen-dently. For any consecutive elements t < t in T i , with t ∈ λ ( v i , ) and t ∈ λ ( v i , ), we add theedge { v i , , v i , } . If p ∞ exists, we add an edge to the graph connecting the vertices correspondingto the last element of T (cid:96) and the first element of the T . Theorem 23 ( PTOPO inside the characteristic box) . Consider a proper parametrization φ of curve C involving polynomials of degree d and bitsize τ , as in Eq.(1), that has no singularities at infinity.Alg. 3 outputs a graph G that, if embedded in R n , is isotopic to C , within the characteristic box B C . It has worst case complexity (cid:101) O B ( d ( n + τ ) + nd τ + n d τ + d ( n τ + n ) + n d τ ) , while its expected complexity is (cid:101) O B ( d τ + nd τ + n d τ + d ( n τ + n ) + n d τ ) . If n = O (1) , then bounds become (cid:101) O B ( N ) , where N = max { d , τ } . Notice that we exclude the parameters of the isolated points. To avoid multiple edges, we make the convention that we add an edge between v i , j , j = ,
2, and an (artificial)intermediate point corresponding to a parameter in ( t , t ). lgorithm 3: PTOPO ( φ ) (Inside the characteristic box) Input:
A proper parametrization φ ∈ Z ( t ) n without singular points at infinity. Output:
Abstract graph G Compute real poles T R P . Compute parameters of ‘special points’ T C , T M , T E , T I . /* Characteristic box */ b ← d ( τ + log d ) , B C ← [ − b , b ] n T B ← parameters that give to intersections of C with B C Construct the set of vertices of G using Lem.22 Sort the list of all the parameters T = [ T C , T M , T E , T B ]. Let T , . . . , T (cid:96) the sublists of T when split at parameters in T R P for every list T i = [ t i , , . . . , t i , k i ] do for j = , . . . , k i − do Add the edge { t i , j , t i , j + } to the graph if p ∞ exists then Add the edge { t , , t (cid:96), k (cid:96) } to the graph Proof.
We count on the fact that φ is continuous in R \ T R P . Thus, for each real interval [ s , t ] with[ s , t ] ∩ T R P = ∅ , there is a parametric arc connecting the points φ ( s ) and φ ( t ). Since for any (sorted)list T i , for i ∈ [ (cid:96) ], the interval defined by the minimum and maximum value of its elements hasempty intersection with T R P , then for any s , t ∈ T i there exists a parametric arc connecting φ ( s )and φ ( t ) and it is entirely contained in B C . If p ∞ exists, then p ∞ is inside B C . Let t , , t (cid:96), k (cid:96) bethe first element of the first list and the last element of the last list. There is a parametric arcconnecting φ ( t , ) with p ∞ and p ∞ with φ ( t (cid:96), k (cid:96) ). So we add the edge { t , , t (cid:96),κ (cid:96) } to G . Then, everyedge of G is embedded to a unique smooth parametric arc and the embedding of G can be triviallycontinuously deformed to C .For the complexity analysis, we know from Lem.18 that steps 1-2 can be performed in wost-case bit complexity (cid:101) O B ( nd + nd τ + d ( n + n τ ) + d ( n τ + n ) + n d τ ) , and in expected bit complexity (cid:101) O B ( d + d ( n + τ ) + d ( n + n τ ) + d ( n τ + n ) + n d τ ) , using a Las Vegas algorithm. From Lemmata 20, 21, and 22 steps 4-5 cost (cid:101) O B ( n ( d τ )).To perform steps 6-7 we must sort all the parameters in T ∪ T R P , i.e., we sort O ( d + nd )algebraic numbers: The parameters that correspond to cusps and extreme points can be expressedas roots of (cid:81) i ∈ [ n ] h i ( t , t ), which is of size ( nd , n τ ). The poles are roots of (cid:81) i ∈ [ n ] q i ( t ), which hassize ( nd , n τ ). The parameters that correspond to multiple points are roots of R s which has size( d , d τ ). At last, parameters in T B are roots of a polynomial of size ( d , d τ ).We can consider all these algebraic numbers together as roots of a single univariate polyno-mial (the product of all the corresponding polynomials). It has degree O ( d + nd ) and bitsize (cid:101) O ( d τ + n τ ). Hence, its separation bound is (cid:101) O ( d τ + nd τ + nd τ + n d τ ). To sort the list of17 igure 1: The left figure is the output of ptopo for the parametric curve ( t + t + t − t − t − , ( t − t + t t − t − t − ) , while the right figureis the output for the curve ( t − t + t − t + t − t + t + t + t + , − t − t − t + t − t − t + tt + t + t + t + ) . Figure 2:
The left figure is the output of ptopo for the parametric curve ( t − t + t − t + , t − t + t − t ) while the right figure is the output for the curve ( t − t + t + t + t − t + , − t − t − t + t − t − t + ) . .all the algebraic numbers we have to perform O ( d + nd ) comparisons and each costs (cid:101) O ( d τ + nd τ + nd τ + n d τ ). Thus, the overall cost for sorting is (cid:101) O B ( d τ + nd τ + n d τ + n d τ + n d τ ).The overall bit complexities in the worst and expected case follow by summing the previousbounds.Following the proof of Thm. 23 we notice that the term d τ in the worst case bound is dueto the introduction of the intersection points of C with B C . For visualizing the curve within B C ,these points are essential and we cannot avoid them. However, if we are interested only in thetopology of C , i.e., the abstract graph G , these points are not important any more. We sketch aprocedure to avoid them and gain a factor of d in the complexity bound:Assume that we have not computed the points on C ∩ B C . We split again the sorted list T = [ T C , T M , T E ] at the real poles, and we add an artificial parameter at the beginning and at theend of each sublist. The rest of the procedure remains unaltered.To verify the correctness of this approach, it su ffi ces to prove that the graph that we obtainby this procedure, is isomorphic to the graph G . It is immediate to see that the latter holds,possibly up to the dissolution of the vertices corresponding to the first and last artificial vertices.Adding these artificial parameters does not a ff ect the overall complexity, since we do not performany algebraic operations. Therefore, the bit complexity of the algorithm is determined by thecomplexity of computing the parameters of the special points (Lem.18), and so we have thefollowing theorem: 18 heorem 24 ( PTOPO and an abstract graph) . Consider a proper parametrization φ of curve C involving polynomials of degree d and bitsize τ , as in Eq. (1), that has no singularities at infinity.Alg. 3 outputs a graph G that, if we embed it in R n , then it is isotopic to C . It has worst casecomplexity (cid:101) O B ( nd + nd τ + d ( n + n τ ) + d ( n τ + n ) + n d τ ) , while its expected complexity is (cid:101) O B ( d + d ( n + τ ) + d ( n + n τ ) + d ( n τ + n ) + n d τ ) , If n = O (1) , then bounds become (cid:101) O B ( N ) , where N = max { d , τ } . Remark 25.
If we are given a box
B ⊂ R n at the input, we slightly modify PTOPO , as follows: Wediscard the parameter values in T C ∪ T M ∪ T E ∪ T I that correspond to points not contained in B . Theset of G’s vertices is constructed similarly. To connect the vertices, we follow the same methodwith a minor modification: For any consecutive elements t < t in a list T i with more than two el-ements, such that t ∈ λ ( v i , ) and t ∈ λ ( v i , ) , we add the edge { v i , , v i , } if and only if φ ( t ) , φ ( t ) are not both on the boundary of B ; or in other words t and t are not both in T B .
6. Isotopic embedding for the special cases of plane and space curves
In this section we elaborate on the isotopic embedding of the output graph G of Thm. 24 forthe case of plane and space parametric curves C . We embed every edge of the abstract graph G to the corresponding parametric arc by sampling many parameter values in the associatedparametric interval and then connect the corresponding points accordingly, in R or R . Thelarger the number of sampled parameters is, the more likely it is for the embedding to be isotopicto C . However, we might need a prohibitive large number of points to sample; their numberis related to the distance between two branches of the curve. We show that by introducing afew additional points, we can replace the parametric arcs of the embedded graph with straightline segments and count on it being isotopic to C . Following closely Alc´azar et al. (2020), if X , Y ⊂ R n are one dimensional, then being isotopic implies that one of them can be deformedinto the other without removing or introducing self-intersections .For plane curves, there is no need to take intermediate points on each parametric arc. Weconsider the embedding of the abstract graph G in R as a straight-line graph ˜ G , i.e., with straightlines for edges, whose vertices are mapped to the corresponding points of the curve. The verticesof ˜ G are all the singular and extreme points with respect to the x - and y - directions. Therefore, theedges of ˜ G correspond to smooth and monotonous parametric arcs and so they cannot intersectbut at their endpoints. The embedding ˜ G is then trivially continuously deformed to C . Thisdiscussion summarizes as follows: Corollary 26 ( PTOPO and isotopic embedding for plane curves) . Consider a proper parametriza-tion φ of a curve C in R involving polynomials of degree d and bitsize τ , as in Eq. (1), that hasno singularities at infinity. Alg. 3 computes an abstract graph whose straight-line embedding in R is isotopic to C in worst case complexity (cid:101) O B ( d + d τ ) . For space curves, a straight-line embedding of G is not guaranteed to be isotopic to C forknots may be present. To overcome this issue we need to segment some edges of G to twoor more edges. To find the extra vertices that we need to add to the graph we follow a commonapproach (Alc´azar and D´ıaz-Toca, 2010; Alc´azar et al., 2020; Kahoui, 2008; Daouda et al., 2008;19 a) (b) Figure 3: (a) Graph of the curve parametrized by φ ( t ) = (cid:0) t − t + t + t − t + t − t + , t + t − t + t − t − t + t − t + , − t − t + t − t + t + t + (cid:1) (green) and its orthogonal projection on the xy-plane (blue). In red, a double point in the projectedcurve with two points in its preimage (apparent singularities). (b) Graph of the curve parametrized by φ ( t ) = (cid:0) − t + t − t − t + − t − t + t − , − t + t − t + t + − t − t + t − , t − t + t − t + t + − t − (cid:1) (green) and its orthogonal projectionon the xy-plane (blue). The space curve does not have a point at infinity, whereas the plane curve has. Cheng et al., 2013) that projects the space curve to a plane one. For a projection defined by themap π : C → R , we write ˜ C = π ( C ). We will ensure in the sequel that the following twoconditions are satisfied:(C1) C has no asymptotes parallel to the direction of the projection.(C2) The map π is birational (Sendra et al., 2008, Def. 2.37).The first condition is to ensure that the point at infinity p ∞ of C exists if and only if the pointat infinity of ˜ C exists and is equal to π ( p ∞ ) (Alc´azar and D´ıaz-Toca, 2010, Lem. 10); see Fig.5afor an instance where this condition is violated. The second condition ensures that only a finitenumber of points on ˜ C have more than one point as a preimage. We call these points apparentsingularities (Daouda et al., 2008); see Fig. 5b. Thus, with this condition we avoid the ”bad”cases where two branches of C project to the same branch of ˜ C . Lemma 27.
Consider a proper parametrization φ of curve C in R involving polynomials ofdegree d and bitsize τ , as in Eq. (1), that has no singularities at infinity. We compute a map π : C → R satisfying conditions (C1), (C2) in worst case complexity (cid:101) O B ( d + d τ ) and using aLas Vegas algorithm in expected complexity (cid:101) O B ( d + d τ ) .Proof. By (Walker, 1978, Thm. 6.5, pg. 146), any space curve can be birationally projected toa plane curve. We choose an integer a uniformly at random from the set { , . . . Kd } , where K = O (1); we explain later in the proof about the size of this set. We define the mapping: π : C → C ( x , y , z ) (cid:55)→ ( x , y + az ) (5)It will be useful in the sequel to regard the application of π ( · ) to φ ( t ) as being performed intwo steps: 20. Change of orthogonal basis of R : Let { , − a , a } be an orthogonal basis of R and A = a − a . In the new basis the curve is parametrized by:( φ ( t ) , φ ( t ) , φ ( t )) · A = ( φ ( t ) , φ ( t ) + a φ ( t ) , − a φ ( t ) + φ ( t ))2. Orthogonal projection onto the first two coordinates: This yields the plane curve parametrizedby ( φ ( t ) , φ ( t ) + a φ ( t )).For a given choice of a we check if conditions (C1) and (C2) are satisfied:For (C1): The direction of the projection is defined by the vector (0 , a , C has noasymptotes parallel to (0 , a ,
1) if and only if the curve with parametrization φ ( t ) · A has no asymp-totes parallel to (0 , a , · A − = (0 , , O ( d ) , (cid:101) O ( τ )) which costs (cid:101) O B ( d + d τ ) worst case.Morover, there are O ( d ) bad values of a for whom (C1) does not hold: for any asymptote of C , there is a unique value of a that maps it to an asymptote parallel to (0 , ,
1) in the new basis.The asymptotes of C are O ( d ) since they occur at the poles of φ ( t ) and at the branches that extendto infinity.For (C2): Since φ ( t ) is proper, π ( · ) is birational if and only if π ( φ ( t )) is also proper. We checkthe properness of this parametrization of ˜ C in (cid:101) O B ( d + d τ ) expected time, using Lem. 5.To find the values of a that result to a ‘bad’ map: Let ˜ h ( s , t ) , ˜ h ( s , t ) the polynomials ofEq. (1) associated to π ( φ ( t )). The parametrization π ( φ ( t )) is proper if and only if gcd( ˜ h ( s , t ) , ˜ h ( s , t )) =
1. If gcd( ˜ h ( s , t ) , ˜ h ( s , t )) (cid:44) R ( s ) = res t ( ˜ h ( s , t ) , ˜ h ( s , t )), we have that R ( s ) = R ( s ) is not always identically zero (e.g., for ˜ h ( s , t ) = t + s , ˜ h ( s , t ) = t + s − R ( s ) = R ( s ) as a polynomial in Z ( a )[ s ]: R ( s ) = c d ( a ) s d + · · · + c ( a ) s + c ( a ) , where c i ∈ Z [ a ] is of size ( d , (cid:101) O ( d τ )), 0 ≤ i ≤ d . The bad values of a ∈ R satisfy then theequation: c d ( a ) + · · · + c ( a ) + c ( a ) = . The polynomial has degree O ( d ) and so there are O ( d ) bad values to avoid. This points to theworst case complexity (cid:101) O B ( d + d τ ).Given a map π computed through Lem. 27, we find the parameters that give the real multiplepoints of ˜ C and its extreme points with respect to the two coordinate axes. We add the corre-sponding vertices to G and we obtain an augmented graph, say G (cid:48) . The straight-line embeddingof G (cid:48) in R is isotopic to ˜ C by Cor. 26, possibly up to the isolated points (Alc´azar and D´ıaz-Toca,2010, Thm. 13 and Lem. 15). Then, by lifting this embedding to the corresponding straight-line graph in R , we obtain an isotopic graph to C (Alc´azar and D´ıaz-Toca, 2010, Thm. 13 andThm. 14). The following theorem summarizes the previous discussion and states its complexity:21 heorem 28 ( PTOPO and isotopic embedding for 3D space curves) . Consider a proper parametriza-tion φ of curve C in R involving polynomials of degree d and bitsize τ , as in Eq. (1), that has nosingularities at infinity. There is an algorithm that computes an abstract graph whose straight-line embedding in R is isotopic to C in worst case complexity (cid:101) O B ( d + d τ ) .Proof. Given a projection map π ( · ) such that (C1) and (C2) hold, correctness follows from theprevious discussion. From Lem. 27, we find such a map in expected complexity (cid:101) O B ( d + d τ ) and (cid:101) O B ( d + d τ ). Using Alg. 2 for π ( φ ( t )) we find the parameters of the extreme and real multiplepoints of ˜ C , say ˜ T , in (cid:101) O B ( d + d τ ) (Lem. 18). Then, we employ Alg. 3 for φ ( t ), by augmentingthe list of parameters that are treated with ˜ T . The last step dominates the complexity and is nota ff ected by the addition of extra parameters to the list since | ˜ T | = O ( d ). The complexity result inThm. 24 allows to conclude. Remark 29.
The complexity in Thm. 28 is no higher than the complexity in Thm. 24, i.e., nocostlier to than the computation of the topology itself.
7. Multiplicities and characterization of singular points
To determine the multiplicity and the character of each real singular point in C we followthe method presented in the series of papers P´erez-D´ıaz (2007); Blasco and P´erez-D´ıaz (2017);Blasco and P´erez-D´ıaz (2019). They provide a complete characterization using resultant compu-tations that applies to curves of any dimension. In the sequel, we present the basic ingredients oftheir approach and we estimate the bit complexity of the algorithm.Let n = H i ( s , t ) = p i ( s ) q i ( t ) − p i ( t ) q i ( s ), for i ∈ [2]. Consider a point p ∈ C given bythe parameter values { s , . . . , s k } , k ≥
1; that is p = φ ( s i ), for all i ∈ [ k ]. The fiber function at p (Blasco and P´erez-D´ıaz, 2019, Def. 2) is F p ( t ) = gcd( H ( s j , t ) , H ( s j , t )) , for any j ∈ [ k ]; it is a univariate polynomial, the real distinct roots of which are the values ofthe parameters that correspond to the point p . When the parametrization of the curve is proper,for any point p on C other than p ∞ it holds that deg( F p ( t )) = mult p ( C ) (Blasco and P´erez-D´ıaz,2019, Cor. 1). Also, when the parameters in { s , . . . , s k } are all real, k equals the number of realbranches that go through p .To classify ordinary and non-ordinary singularities we proceed as follows: We say that asingularity is ordinary if it is at the intersection of smooth branches only and the tangents toall branches are distinct (Walker, 1978, p.54). In all the other cases, we call it a non-ordinary singularity. Non-ordinary singularities may have other singularities in their neighborhood. For apoint p ∈ C the delta invariant δ p is a nonnegative integer that measures the number of doublepoints concentrated around p . We can compute it by taking into account the multiplicities of p and its neighboring singularities. Blasco and P´erez-D´ıaz (2019) consider three di ff erent types ofnon-ordinary singularities and use the delta invariant to distinguish them. In particular:1. If k = mult p ( C ) and 2 δ p = mult p ( C )(mult p ( C ) − p is an ordinary singularity.2. If k < mult p ( C ) and 2 δ p = mult p ( C )(mult p ( C ) − p is a type I non-ordinary singu-larity. 22. If k = mult p ( C ) and 2 δ p > mult p ( C )(mult p ( C ) − p is a type II non-ordinarysingularity.4. If k < mult p ( C ) and 2 δ p > mult p ( C )(mult p ( C ) − p is a type III non-ordinarysingularity.Blasco and P´erez-D´ıaz (2019) compute the delta invariant, δ p , using the formula δ p = k (cid:88) j = (cid:88) t s.t. h ( s j , t ) = h ( s j , t ) = Int h , h ( s j , t ) , (6)where Int h , h ( α, β ) is the intersection multiplicity of the coprime polynomials h ( s , t ) and h ( s , t )at a point ( α, β ). Using a well known result, e.g., (Fulton, 1984, 1.6) as stated in (Bus´e et al.,2005, Prop. 5), we can compute the intersection multiplicities using resultant computations. Let R ( s ) = res t ( h ( s , t ) , h ( s , t )). For a root α of R ( s ), its multiplicity µ ( a ) is equal to the sum of theintersection multiplicities of solutions of the system of { h , h } in the form ( α, t ), and so µ ( α ) = (cid:88) t s.t. h ( α, t ) = h ( α, t ) = Int h , h ( α, t ) . (7)Therefore, from Eq. (6) and Eq. (7), we conclude that: δ p = k (cid:88) j = µ ( s j ) . (8)For space curves, we birationally project C to a plane curve (Walker, 1978, Thm. 6.5, pg. 146).The multiplicities of the singular points of C and their delta invariant are preserved under the bi-rational map realizing the projection (Blasco and P´erez-D´ıaz, 2019, Prop. 1 and Cor. 4). Thepseudo code of the algorithm appears in Alg. 4. Theorem 30.
Let C be a curve with a proper parametrization φ ( t ) as in Eq. (1), that has nosingularities at infinity. Alg. 4 computes the singular points of C , their multiplicity and character(ordinary / non-ordinary) in (cid:101) O B ( d + d τ ) worst-case complexity when n = and for n > in expected complexity (cid:101) O B ( d + d ( n + τ ) + d ( n + n τ ) + d ( n τ + n ) + n d τ ) . Proof.
We compute the parameters that give the singular points of C using Alg. 2 in (cid:101) O B ( d + d ( n + τ ) + d ( n + n τ ) + d ( n τ + n ) + n d τ ) when n >
2, which becomes worst-case when n = n >
2, lines 5-8 compute a birational projection of C to a plane curve parametrized by˜ φ ( t ) (note that the projection is birational is and only if ˜ φ ( t ) is proper since φ ( t ) is proper). Theexpected complexity of this is (cid:101) O B ( d + nd τ ) by slightly adapting the proof of Lem. 27.We can group the parameter values that give singular points using Lem. 22 in O ( d + nd )arithmetic operations. We have that gcd( H ( s , t ) , H ( s , t )) = s − t . For h ( s , t ) = H ( s , t ) / ( s − t ), h ( s , t ) = H ( s , t ) / ( s − t ), we compute a triangular decomposition of the system { h ( s , t ) = h ( s , t ) = } ; it consists of the systems { ( A i ( s ) , B i ( s , t )) } i ∈I . For any root α of A i , B i ( α ; t ) is of23 lgorithm 4: CharacterizeSingularPoints ( φ, S ) Input:
Proper parametrization φ ∈ Z ( t ) n without singularity at infinity, as in Eq. (1) Output:
Multiplicities and characterization of points S ←
Special Points ( φ ) if n = then Compute polynomials H ( s , t ) , H ( s , t ) for φ else repeat Choose integers a , . . . , a n at random from { , . . . , Kd n } , where K = O (1). ˜ φ ( t ) ← ( φ ( t ) , φ ( t ) + a φ ( t ) + · · · + a n φ n ( t )) until ˜ φ ( t ) is proper ; Compute polynomials H ( s , t ) , H ( s , t ) for ˜ φ for p ∈ S do M p ← { s ∈ R : φ ( s ) = p } // parameters that give the same point p k ← | M p | // number of real branches that go through p Take s ∈ M p14 mult p ( C ) ← deg(gcd( H ( s , t ) , H ( s , t ))) // multiplicity// compute delta invariant δ p ← for s j ∈ M p do µ ( s j ) ← multiplicity of s j as a root of res t ( H ( s , t ) / ( s − t ) , H ( s , t ) / ( s − t )) δ p ← δ p + µ ( s j ) δ p ← δ p / if k = mult p ( C ) and δ p = mult p ( C )( mult p ( C ) − then return p is an ordinary singularity else if k < mult p ( C ) and δ p = mult p ( C )( mult p ( C ) − then return p is a type I non-ordinary singularity else if k = mult p ( C ) and δ p > mult p ( C )( mult p ( C ) − then return p is a type II non-ordinary singularity else if k < mult p ( C ) and δ p > mult p ( C )( mult p ( C ) − then return p is a type III non-ordinary singularitydegree i and equals gcd( h ( α ; t ) , h ( α ; t )) (up to a constant factor). By (Bouzidi et al., 2015a,Prop. 16), the triangular decomposition is computed in (cid:101) O B ( d + d τ ) worst case .For a singular point p and s ∈ M p : To compute the degree of gcd( H ( s ) , H ( s )), since s is a root of res t ( h ( s , t ) , h ( s , t )), it su ffi ces to find for which i ∈ I it holds A i ( s ) = s is given in isolated interval representation as a root of res t ( h ( s , t ) , h ( s , t )). However, we can isolate the roots of all A i by taking care of the isolat-ing invervals to become small enough such that the isolating interval of s intersects only one In a Las Vegas setting this computation could be reduced to (cid:101) O B ( d + d τ ), but since it does not a ff ect the totalcomplexity we chose not to expand onto this. a) (b) Figure 4:
Output of
PTOPO for (a) a Lissajous curve parametrized by (cid:0) ( − t − t + , − t ( t − t − t + t + t + t + t + , − (16( t − t + t − t ( t − t + t − t + t + t + t + t + t + t + t + t + (cid:1) . The multiple points are indicated by the redsquares. (b) a space curve parametrized by (cid:0) − t t + , − t (3 t − t + t − t − t + t + t − t + t + , t (cid:1) . The apparent singularities are indicated bythe black squares. of them. Because they are roots of the same polynomial this cannot exceed the complexity ofisolating the roots of res t ( h ( s , t ) , h ( s , t )). This is can be done in (cid:101) O B ( d + d τ ). In the same bitcomplexity we have the multiplicity µ ( s j ) of s j as a root of res t ( h , h ).
8. Implementation and examples
PTOPO is implemented in maple . The current public version works for plane curves and wewill release the version for the 3D curves in the near future. A typical output appears in Figures 1,2, 4, and 5. Besides the visualization the software computes all the points of interest of curve(singular, extreme, etc) in isolating interval representation as well as in suitable floating pointapproximations.We build upon the real root isolation routines of maple ’s RootFinding library and the slv package (Diochnos et al., 2009), to use a certified implementation of general purpose exact com-putations with one and two real algebraic numbers, like comparison and sign evaluations, as wellas exact (bivariate) polynomial solving.
PTOPO computes the topology and visualizes parametric curves in two (and in the near futurein three) dimensions. For a given parametric representation of a curve,
PTOPO computes thespecial points on the curve, the characteristic box, the corresponding graph, and then it visualizesthe curve (inside the box). The computation, in all examples from literature we tested, takes lessthan a second in a MacBook laptop, running maple https://webusers.imj-prg.fr/~christina.katsamaki/ptopo/ This project has received funding from the European Union’s Horizon 2020 research and innovation program under theMarie Skłodowska-Curie grant agreement No 754362. a) (b) Figure 5:
Output of
PTOPO for a curve parametrized by (a) (cid:0) − t + t − t − t + − t − t + t − , − t ∗ t − t − t − − t − t + t − , − t ∗ t − t − t − − t − t + t − (cid:1) (b) (cid:0) − t + t + , ( − t + t ( t + , ( − t + t ( t + (cid:1) . Examples are taken from Alc´azar and D´ıaz-Toca (2010). Acknowledgments.
We acknowledge partial support by the Fondation Math´ematique JacquesHadamard through the PGMO grand ALMA, by ANR JCJC GALOP (ANR-17-CE40-0009), bythe PHC GRAPE, by the projects 118F321 under the program 2509, 118C240 under the program2232, and 117F100 under the program 3501 of the Scientific and Technological Research Councilof Turkey.
References
Abhyankar, S. S., Bajaj, C. J., 1989. Automatic parameterization of rational curves and surfaces IV: Algebraic spacecurves. ACM Trans. Graph. 8 (4), 325–334.URL https://doi.org/10.1145/77269.77273
Alberti, L., Mourrain, B., Wintz, J., 2008. Topology and arrangement computation of semi-algebraic planar curves.CAGD 25 (8), 631 – 651.URL
Alc´azar, J. G., Caravantes, J., D´ıaz, G. M., Tsigaridas, E., 2020. Computing the topology of a plane or space hyperellipticcurve. Computer Aided Geometric Design 78, 101830.URL
Alc´azar, J. G., D´ıaz-Toca, G. M., 2010. Topology of 2D and 3D rational curves. CAGD 27 (7), 483 – 502.URL
Basu, S., Pollack, R., M-F.Roy, 2003. Algorithms in Real Algebraic Geometry. Vol. 10 of Algorithms and Computationin Mathematics. Springer-Verlag.Bernardi, A., Gimigliano, A., Id`a, M., 09 2016. Singularities of plane rational curves via projections. J. Symb. Comput.Blasco, A., P´erez-D´ıaz, S., 2017. Resultants and singularities of parametric curves. arXiv.URL https://arxiv.org/abs/1706.08430s
Blasco, A., P´erez-D´ıaz, S., 2019. An in depth analysis, via resultants, of the singularities of a parametric curve. CAGD68, 22–47.Boissonnat, J.-D., Teillaud, M. (Eds.), 2006. E ff ective Computational Geometry for Curves and Surfaces. Springer-Verlag, Mathematics and Visualization.Bouzidi, Y., Lazard, S., Moroz, G., Pouget, M., Rouillier, F., Sagralo ff , M., 02 2015a. Improved algorithms for solvingbivariate systems via rational univariate representations.Bouzidi, Y., Lazard, S., Moroz, G., Pouget, M., Rouillier, F., Sagralo ff , M., 2016. Solving bivariate systems usingRational Univariate Representations. J. Complexity 37, 34–75.URL https://hal.inria.fr/hal-01342211 ouzidi, Y., Lazard, S., Pouget, M., Rouillier, F., 2013. Rational univariate representations of bivariate systems andapplications. In: Proc. 38th Int’l Symp. on Symbolic and Algebraic Computation. ISSAC ’13. ACM, NY, USA, pp.109–116.URL http://doi.acm.org/10.1145/2465506.2465519 Bouzidi, Y., Lazard, S., Pouget, M., Rouillier, F., 2015b. Separating linear forms and rational univariate representationsof bivariate systems. J. Symb. Comput., 84–119.URL https://hal.inria.fr/hal-00977671
Bus´e, L., 2014. Implicit matrix representations of rational B´ezier curves and surfaces. Computer-Aided Design 46, 14–24.URL https://hal.inria.fr/hal-00847802
Bus´e, L., D’Andrea, C., 2012. Singular factors of rational plane curves. J. Algebra 357, 322–346.Bus´e, L., Khalil, H., Mourrain, B., 2005. Resultant-based methods for plane curves intersection problems. In: Ganzha,V. G., Mayr, E. W., Vorozhtsov, E. V. (Eds.), Computer Algebra in Scientific Computing. Springer Berlin Heidelberg,Berlin, Heidelberg, pp. 75–92.Bus´e, L., Laroche, C., Yıldırım, F., 2019. Implicitizing rational curves by the method of moving quadrics. Computer-Aided Design 114, 101–111.Bus´e, L., Luu Ba, T., Oct. 2010. Matrix-based Implicit Representations of Rational Algebraic Curves and Applications.CAGD 27 (9), 681–699.URL https://hal.inria.fr/inria-00468964
Caravantes, J., Fioravanti, M., Gonzalez-Vega, L., Necula, I., 2014. Computing the topology of an arrangement ofimplicit and parametric curves given by values. In: Gerdt, V. P., Koepf, W., Seiler, W. M., Vorozhtsov, E. V. (Eds.),Computer Algebra in Scientific Computing. Springer, Cham, pp. 59–73.Cheng, J.-S., Jin, K., Lazard, D., 2013. Certified rational parametric approximation of real algebraic space curves withlocal generic position method. Journal of Symbolic Computation 58, 18 – 40.URL
Chionh, E.-W., Sederberg, T. W., 2001. On the minors of the implicitization b´ezout matrix for a rational plane curve.CAGD 18 (1), 21 – 36.Cox, D., Kustin, A., Polini, C., Ulrich, B., 02 2011. A study of singularities on rational curves via syzygies. Memoirs ofthe American Mathematical Society 222.Daouda, D. N., Mourrain, B., Ruatta, O., 2008. On the computation of the topology of a non-reduced implicit space curve.In: Proceedings of the Twenty-First International Symposium on Symbolic and Algebraic Computation. ISSAC ’08.Association for Computing Machinery, New York, NY, USA, p. 47–54.URL https://doi.org/10.1145/1390768.1390778
Diatta, D. N., Diatta, S., Rouillier, F., Roy, M.-F., Sagralo ff , M., 2018. Bounds for polynomials on algebraic numbersand application to curve topology. arXiv preprint arXiv:1807.10622.Diochnos, D. I., Emiris, I. Z., Tsigaridas, E. P., 2009. On the asymptotic and practical complexity of solving bivariatesystems over the reals. J. Symb. Comput. 44 (7), 818–835, (Special issue on ISSAC 2007).Farouki, R. T., Giannelli, C., Sestini, A., 2010. Geometric design using space curves with rational rotation-minimizingframes. In: Dæhlen, M., Floater, M., Lyche, T., Merrien, J.-L., Mørken, K., Schumaker, L. L. (Eds.), MathematicalMethods for Curves and Surfaces. Springer, pp. 194–208.Fulton, W., 1969. Algebraic Curves. An Introduction to Algebraic Geometry. Addison Wesley.Fulton, W., 1984. Introduction to Intersection Theory in Algebraic Geometry. Cbms Regional Conference Series inMathematics. American Mathematical Society.Gao, X.-S., Chou, S.-C., 1992. Implicitization of rational parametric equations. J. Symb. Comput. 14 (5), 459 – 470.URL Gutierrez, J., Rubio, R., Sevilla, D., 2002a. On multivariate rational function decomposition. J. Symb. Comput. 33 (5),545 – 562.URL
Gutierrez, J., Rubio, R., Yu, J.-T., 08 2002b. D-resultant for rational functions. Proc. American Mathematical Society130.Jia, X., Shi, X., Chen, F., 2018. Survey on the theory and applications of µ -bases for rational curves and surfaces. J.Comput. Appl. Math. 329, 2–23.Kahoui, M. E., 2008. Topology of real algebraic space curves. J. Symb. Comput. 43 (4), 235 – 258.URL Katsamaki, C., Rouillier, F., Tsigaridas, E., Zafeirakopoulos, Z., 2020a. On the geometry and the topology of parametriccurves. In: Proc. 45th Int’l Symposium on Symbolic and Algebraic Computation. ISSAC ’20. ACM, NY, USA, p.281–288.Katsamaki, C., Rouillier, F., Tsigaridas, E. P., Zafeirakopoulos, Z., 2020b. PTOPO: a maple package for the topology ofparametric curves. ACM Commun. Comput. Algebra 54 (2), 49–52. obel, A., Sagralo ff , M., 08 2014. On the complexity of computing with planar algebraic curves. J. Complexity 31.Li, Y.-M., Cripps, R. J., 1997. Identification of inflection points and cusps on rational curves. CAGD 14 (5), 491 – 497.URL Lickteig, T., Roy, M.-F., Mar. 2001. Sylvester–Habicht sequences and fast Cauchy index computation. J. Symb. Comput.31 (3), 315–341.URL https://doi.org/10.1006/jsco.2000.0427
Mahler, K., 1962. On some inequalities for polynomials in several variables. J. London Mathematical Society 1 (1),341–344.Manocha, D., Canny, J. F., 1992. Detecting cusps and inflection points in curves. CAGD 9 (1), 1 – 24.URL
Pan, V., 2002. Univariate polynomials: Nearly optimal algorithms for numerical factorization and rootfinding. J. Symb.Comput. 33 (5), 701–733.Pan, V., Tsigaridas, E., 2017. Accelerated approximation of the complex roots and factors of a univariate polynomial.Theor. Computer Science 681, 138 – 145.URL
Park, H., 08 2002. E ff ective computation of singularities of parametric a ffi ne curves. J. Pure and Applied Algebra 173,49–58.P´erez-D´ıaz, S., 2006. On the problem of proper reparametrization for rational curves and surfaces. CAGD 23 (4), 307–323.P´erez-D´ıaz, S., 2007. Computation of the singularities of parametric plane curves. J. Symb. Comput. 42 (8), 835 – 857.URL Recio, C. A. T., 02 2007. Plotting missing points and branches of real parametric curves. Applicable Algebra in Engi-neering, Communication and Computing 18.Rubio, R., Serradilla, J., V´elez, M., 2009. Detecting real singularities of a space curve from a real rational parametriza-tion. J. Symb. Comput. 44 (5), 490 – 498.URL
Sch¨onhage, A., 1988. Probabilistic computation of integer polynomial gcds. J. Algorithms 9 (3), 365 – 371.URL
Sederberg, T. W., May 1986. Improperly parametrized rational curves. CAGD 3 (1), 67–75.URL
Sederberg, T. W., Chen, F., 1995. Implicitization using moving curves and surfaces. In: Proc. of the 22nd AnnualConference on Computer Graphics and Interactive Techniques. SIGGRAPH ’95. NY, USA, pp. 301–308.URL https://doi.org/10.1145/218380.218460
Sendra, J. R., 2002. Normal parametrizations of algebraic plane curves. J. Symb. Comput. 33, 863–885.Sendra, J. R., Winkler, F., Apr. 1999. Algorithms for rational real algebraic curves. Fundam. Inf. 39 (1,2), 211–228.URL http://dl.acm.org/citation.cfm?id=2378083.2378093
Sendra, J. R., Winkler, F., P´erez-D´ıaz, S., 2008. Rational algebraic curves. Algorithms and Computation in Mathematics22.Strzebonski, A., Tsigaridas, E., 2019. Univariate real root isolation in an extension field and applications. J. Symb.Comput. 92, 31 – 51.URL van den Essen, A., YU, J.-T., 01 1997. The D-resultant, singularities and the degree of unfaithfulness. Proc. AmericanMathematical Society 125.von zur Gathen, J., Gerhard, J., 2013. Modern computer algebra, 3rd Edition. Cambridge University Press.Walker, R. J., 1978. Algebraic curves. Springer-Verlag.van den Essen, A., YU, J.-T., 01 1997. The D-resultant, singularities and the degree of unfaithfulness. Proc. AmericanMathematical Society 125.von zur Gathen, J., Gerhard, J., 2013. Modern computer algebra, 3rd Edition. Cambridge University Press.Walker, R. J., 1978. Algebraic curves. Springer-Verlag.