Efficient topology determination of implicitly defined algebraic plane curves

Computer Aided Geometric Design 19 (2002) 719–743 www.elsevier.com/locate/cagd

Efficient topology determination of implicitly defined algebraic plane curves Laureano Gonzalez-Vega ∗,1 , Ioana Necula 1 Departamento de Matematicas, Estadistica y Computacion Facultad de Ciencias, Universidad de Cantabria, Santander 39005, Spain Received 25 August 2001

Abstract This paper is devoted to present a new algorithm computing in a very efficient way the topology of a real algebraic plane curve defined implicitly. This algorithm proceeds in a seminumerical way by performing a symbolic preprocessing which allows later to accomplish the numerical computations in a very accurate way.  2002 Elsevier Science B.V. All rights reserved. Keywords: Algebraic curves; Seminumerical algorithms; Topology computation

1. Introduction Many important problems in Computer Aided Geometric Design are reduced to the computation of the graph of an algebraic plane curve presented implicitly. For example if we want to section the surface Y (s, t) Z(s, t) X(s, t) , y= , z= , s, t ∈ [0, 1], x= W (s, t) W (s, t) W (s, t) (defined by the real polynomials X(s, t), Y (s, t), Z(s, t) and W (s, t)) with respect to the plane x = x0 then we have two possibilities: either we draw into the unit square [0, 1] × [0, 1] the algebraic plane curve defined by X(s, t) − x0 W (s, t) = 0 and then this drawing is lifted to the considered surface or, if the implicit equation H (x, y, z) of the considered surface is available, the lifting procedure can be avoided by merely computing the graph * Corresponding author.

E-mail addresses: [email protected] (L. Gonzalez-Vega), [email protected] (I. Necula). 1 Partially supported by DGES PB98-0713-C02-02 and by FEDER project 1FD-97-0409.

0167-8396/02/$ – see front matter  2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 3 9 6 ( 0 2 ) 0 0 1 6 7 - X

720

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

Fig. 1. 2y 3 − (3x − 3)y 2 − (3x 2 − 3x)y − x 3 = 0.

of the algebraic plane curve H (x0 , y, z) = 0. But both possibilities require to draw an algebraic curve implicitly defined. This paper is devoted to present a new seminumerical algorithm, and its implementation into the Computer Algebra System Maple, computing, for a given bivariate polynomial f (x, y) in Z[x, y], a linear graph representing the real algebraic plane curve defined by the considered polynomial even if the considered curve has complicated singularities. The main purpose of this algorithm is to compute in the fastest way correct qualitative information about the set f (x, y) = 0 that could be of further use when drawing the considered curve. For example, for the curve C defined by the polynomial
 f (x, y) = 2y 3 − (3x − 3)y 2 − 3x 2 − 3x y − x 3 , the real drawing of C appears to the left in Fig. 1 and the corresponding graph appears to the right of the same figure. A more complicated case appears when considering the polynomial g(x, y) = 279756x + 279936xy 4 − 559692y 2 x − 15583x 3 + 217x 5 +

130218 2 23039 4 x − x 5 10

359 6 3726432 2 12947 2 4 x + 370656y 4 − y − 72774y 2 x 2 + y x + 1296y 6 10 5 5 37333439 + 46728y 4 x 2 + 15588y 2 x 3 + 100 since when it is drawn (Fig. 2, left), it is difficult to conclude what is its real topological shape (see Fig. 2, right). The problem of computing the graph (even topologically) of an algebraic plane curve defined implicitly has received a special attention, first, from Computer Aided Geometric Design since it is a basic subproblem appearing very often in practice (see, for example, (Bajaj et al., 1988; Keyser et al., 2000, 1999a, 1999b)) and, second, from Computer Algebra since it has been responsible of many advances regarding subresultants, symbolic real root counting, infinitesimal computations, etc. From the seminal papers (Arnon and McCallum, 1988; Arnborg and Feng, 1988; Gianni and Traverso, 1983; Roy, 1996; Sakkalis, 1991), the interested reader can see in (Cellini et al., 1991; Cucker et al., 1991; Roy and +

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

721

Fig. 2. g(x, y) = 0.

Szpirglas, 1990; Feng, 1992; Gonzalez-Vega and El Kahoui, 1996; Hong, 1996) how the theoretical and practical complexities of the algorithms dealing with this problem have been dramatically improved. The usual strategy to compute the graph (even topologically) of a real algebraic plane curve defined implicitly by a polynomial f (x, y) ∈ R[x, y] proceeds in the following way: • Step I: Computation of the discriminant of f with respect to y, Φ0 (x), and determination of its real roots, α1 < · · · < αr . • Step II: For every αi , computation of the real roots of f (αi , y), βi,1 < · · · < βi,si . • Step III: For every αi and βi,j computation of the number of half–branches to the right and to the left of the point (αi , βi,j ). All this information is enough to characterize the required topological information of the curve defined by f (x, y). The proposed method avoids the numerical problems (as clearly quoted in (Hong, 1996)) arising from the computation of the roots of every f (αi , y), which has always multiple roots, by performing a generic linear change of variables, if it is needed, in order to have the following condition verified for every α ∈ R:   ∂f # β ∈ C: f (α, β) = 0, (α, β) = 0
1. ∂y This assures that for every αi real root of Φ0 (x), there is at most one critical point of the curve (a singular point or a point with a vertical tangent) in the vertical line x = αi and that the y-coordinate of such a point can be rationally described in terms of αi (as it will be shown in Section 3). Moreover this allows to symbolically construct, from every f (αi , y), a squarefree polynomial g(αi , y) whose real roots are to be computed in order to finish with the so called Step II. Step III is thus accomplished by merely computing the number of real roots of the squarefree polynomials f (γi , y) (i ∈ {0, 1, . . . , r}) with γ0 = −∞, γr = ∞ and γi any real number in the open interval (αi , αi+1 ) which motivates the squarefreness of f (γi , y).

722

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

These computations provide a graph of the considered curve which is very helpful when this curve is going to be traced numerically since we know exactly how to proceed when close to a complicated point. The performed experimentation shows that the proposed method is faster than the existing ones and provides a first example where breaking the sparness of the initial polynomial equation in order to get a simpler geometry makes the computations to be performed easier. The paper is divided in five sections. In the second section we recall definitions and main properties of Sturm–Habicht sequence related with the real root counting problem. Third section is devoted to present the definition of curve in generic position together with the main properties of the critical points for curves in generic position. Fourth section presents a detailed description of the algorithm together with the complete description of its application to a concrete example. Fifth section shows how the implementation in the Computer Algebra System Maple was accomplished together with a description of the performed experimentation. Finally last section drafts several conclusions together with several future extensions of the algorithm presented in the fourth section.

2. Algebraic preliminaries: Sturm–Habicht sequence This section is devoted to introduce the definition of the Sturm–Habicht sequence and its main properties related with the real root counting problem. Sturm–Habicht sequence was introduced in (Gonzalez-Vega et al., 1990, 1994); proofs of the results summarized into this section can be found in (Gonzalez-Vega et al., 1990, 1994) or (Gonzalez-Vega et al., 1998). Definition 2.1. Let P , Q be polynomials in R[x] with deg(P ) = p and deg(Q) = q: P=

p 

Q=

k

ak x ,

k=0

q 

bk x k .

k=0

If i ∈ {0, . . . , inf(p, q) − 1} we define the polynomial subresultant associated to P and Q of index i as follows: Sresi (P , Q) =

i 

Mji (P , Q)x j ,

j =0

where every Mji (P , Q) is the determinant of the matrix built with the columns 1, 2, . . . , p + q − 2i − 1 and p + q − i − j in the matrix:  a

p

   mi (P , Q) =   bq  

p+q−i

... .. . ... .. .

 a0

..

ap b0

. ...

    a0      

q −i .

p−i . bq . . . b0 i The determinant Mi (P , Q) will be called ith principal subresultant coefficient and denoted by sresi (P , Q). ..

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

723

Next definition introduces Sturm–Habicht sequence associated to P as the subresultant sequence for P and P  modulo certain sign changes. Definition 2.2. Let P be a polynomial in R[x] with p = deg(P ). If we write δk = (−1)k(k+1)/2 for every integer k, the Sturm–Habicht sequence associated to P is defined as the list of polynomials {StHaj (P )}j =0,...,p where StHap (P ) = P , StHap−1 (P ) = P  and for every j ∈ {0, . . . , p − 2}: StHaj (P ) = δp−j −1 Sresj (P , P  ). For every j in {0, . . . , p} the principal j th Sturm–Habicht coefficient is defined as: 
sthaj (P ) = coefj StHaj (P ) , i.e., the coefficient of x j in the polynomial StHaj (P ). It is important to quote that the discriminant of P is equal to the polynomial stha0 (P ) modulo the leading coefficient of P : stha0 (P ) = ap discriminant(P ). Sign counting on the principal Sturm–Habicht coefficients provides a very useful information about the real roots of the considered polynomial. Next definitions show which are the sign counting functions to be used in the sequel (see (Gonzalez-Vega et al., 1990, 1994) or (Gonzalez-Vega et al., 1998)). Definition 2.3. Let I = {a0 , a1 , . . . , an } be a list of nonzero elements in R. We define: • V(I) as the number of sign variations in the list {a0 , a1 , . . . , an }, • P(I) as the number of sign permanences in the list {a0 , a1 , . . . , an }. Definition 2.4. Let a0 , a1 , . . . , an be elements in R with a0 = 0 and with the following distribution of zeros: I = {a0 , a1 , . . . , an } k

k

 1  2  = a0 , . . . , ai1 , 0, . . . , 0, ai1 +k1 +1 , . . . , ai2 , 0, . . . , 0, ai2 +k2 +1 , , ai3 , 0, . . . . . . . . . , 0, kt

   ait−1 +kt−1 +1 , . . . , ait , 0, . . . , 0 , where all the ai ’s that have been written are not 0. We define i0 + k0 + 1 = 0 and C(I) =

t t −1  
 

εis , P {ais−1 +ks−1 +1 , . . . , ais } − V {ais−1 +ks−1 +1 , . . . , ais } + s=1

where εis =



s=1

0 if ks is odd,
ais +ks +1  ks /2 if ks is even. (−1) sign ai s

724

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

Next the relation between the real zeros of a polynomial P ∈ R[x] and the polynomials in the Sturm– Habicht sequence of P is presented. Its proof can be found in (Gonzalez-Vega et al., 1990, 1994) or (Gonzalez-Vega et al., 1998). Theorem 2.5. If P is a polynomial in R[x] with p = deg(P ) then:   
 C sthap (P ), . . . , stha0 (P ) = # α ∈ R: P (α) = 0 . In particular we have that the number of real roots of P is determined exactly by the signs of the last p − 1 determinants sthai (P ) (the first two ones are lcof(P ) and p lcof(P ) with lcof(P ) denoting the leading coefficient of P ). Moreover, as these determinants are the formal leading coefficients of the subresultant sequence for P and P  (modulo some sign changes), the greatest common divisor of P and P  is obtained as a by-product together with the following equivalence:  stha0 (P ) = · · · = sthai−1 (P ) = 0  (1) StHai (P ) = gcd(P , P ) ⇔ sthai (P ) = 0 which will be very useful in Section 3. Note that when stha0 (P ) = 0, the greatest common divisor of P and P  is trivial: in other words P is squarefree (i.e., without multiple roots). The definition of Sturm–Habicht sequence through determinants allows to perform computations dealing with real roots in a generic way: if P and Q are two polynomials with parametric coefficients whose degrees do not change after specialization we can compute the Sturm–Habicht sequence for P and Q without specializing the parameters and the result is always good after specialization (modulo the condition over the degrees). In the case considered here, this is specially useful since the Sturm–Habicht sequence of f (x, y) will be computed considering y as the main variable and then, by merely replacing x by α, the Sturm–Habicht sequence of f (α, y) will be obtained. This behaviour is studied carefully in (Gonzalez-Vega et al., 1990, 1994) where Sturm–Habicht sequence definition is modified in order to get good specialization properties even in cases where the degrees change after specialization. This is not true when using Sturm sequences (the computation of the remainders makes to appear denominators which can vanish after specialization) or negative polynomial remainder sequences (with fixed degree for P the sequence has not always the same number of elements); see (Loos, 1982; Gonzalez-Vega et al., 1990, 1994) for a more detailed explanation. To quote finally, with respect to Sturm–Habicht sequences, that when the coefficients of P belong to Z or Z[x], the most efficient way of computing the sthaj (P )’s is through Subresultant Algorithm by means of O(deg(P )2 ) arithmetic operations in Z or Z[x]. The size of the involved integers and degrees of intermediate results are well bounded (see (Gonzalez-Vega et al., 1990, 1994), and the more recent (Lombardi et al., 2000)). When the coefficients of P are floating point real numbers the Subresultant Algorithm is no longer applicable and the Definition 2.2 should be directly applied. 3. Critical and regular points for curves in generic position The notion of curve in generic position is introduced and shown to have a very good computational behaviour when considering critical points. We show that the critical point set for these curves is very easy to describe by means of the Sturm–Habicht sequence and that it is possible to get a change of coordinates such that the curve, with respect to the new coordinates, is in generic position. Next the definitions of critical, singular and regular point are introduced.

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

725

Definition 3.1. Let f (x, y) ∈ R[x, y] and C(f ) be the real algebraic plane curve defined by f :   C(f ) = (α, β) ∈ R2 : f (α, β) = 0 . A point (α, β) ∈ C2 is called: • a critical point of C(f ) if ∂f (α, β) = 0; f (α, β) = 0, ∂y • a singular point of C(f ) if ∂f ∂f (α, β) = 0, (α, β) = 0; f (α, β) = 0, ∂y ∂x • a regular point of C(f ) if f (α, β) = 0 and it is not critical. Next two propositions show how the properties of Sturm–Habicht sequences can be used to deal with the set of critical points of a real algebraic plane curve. It is well known that the real roots of stha0 (f ) (if considering f (x, y) as a polynomial in y), which are the real roots of the discriminant of f with respect to y, contain the x-coordinates of the critical points of f . Note that if x = α is a vertical asymptote of C(f ) then α is a real root of stha0 (f ) and that the x-coordinates of complex, non-real critical points of C(f ) with x-coordinate real are also real roots of stha0 (f ). Proposition 3.2. Let f (x, y) ∈ R[x, y] be a squarefree polynomial such that lcofy (f ) ∈ R (or without real roots when belonging to R[x] − R) and α ∈ R. If hj (x) = sthaj (f ) and Hj (x, y) = StHaj (f ) (computed with respect to y) then • Critical point control: h0 (α) = 0, . . . , hk−1 (α) = 0, hk (α) = 0

⇒

  ∂f gcd f (α, y), (α, y) = Hk (α, y). ∂y

• Number of points in a vertical line:   
 # β ∈ R: f (α, β) = 0 = C hd (α), . . . , h0 (α) . Proof. The implication providing the critical point control is directly derived from equivalence (1) in Section 2. The formula for the number of points in a vertical line is a direct consequence of Theorem 2.5. ✷ Note that the polynomials hj (x) defined in the previous proposition are exactly, modulo the sign, the principal subresultants of f and its partial derivative with respect to y:   ∂f . hj (x) = δp−j −1 sresj f, ∂y They are going to play a fundamental role in what follows. First, next proposition presents a nice description of the set of the critical points for the real algebraic plane curve defined by the polynomial f (x, y).

726

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

Proposition 3.3. Let f (x, y) ∈ R[x, y] be a squarefree polynomial of degree d in y such that lcofy (f ) ∈ R (or without real roots when belonging to R[x] − R). If hj (x) = sthaj (f ) and Hj (x, y) = StHaj (f ) (computed with respect to y) then  d−1    ∂f 2 (α, β) ∈ C2 : Γk (α) = 0, Hk (α, β) = 0 , (α, β) ∈ C : f (α, β) = 0, (α, β) = 0 = ∂y k=1 where the polynomials Γk (x) are defined inductively by the following relations: h0 (x) , gcd(h0 (x), h0 (x))
 Φ1 (x) = gcd Φ0 (x), h1 (x) ,

Φ0 (x) =


 Φ2 (x) = gcd Φ1 (x), h2 (x) ,
 Φ3 (x) = gcd Φ2 (x), h3 (x) , .. .


 Φd−2 (x) = gcd Φd−3 (x), hd−2 (x) ,
 Φd−1 (x) = gcd Φd−2 (x), hd−1 (x) ,

Φ0 (x) , Φ1 (x) Φ1 (x) , Γ2 (x) = Φ2 (x) Φ2 (x) , Γ3 (x) = Φ3 (x) .. . Φd−3 (x) , Γd−2(x) = Φd−2 (x) Φd−2 (x) . Γd−1(x) = Φd−1 (x) Γ1 (x) =

Proof. The description of the set of critical points to be proved comes from the different possibilities for (α, y) when α is a root of the discriminant. the degree of the greatest common divisor of f (α, y) and ∂f ∂y Note that, due to equivalence (1) in Section 2, this degree is determined by the vanishing of the different hj (x) since the roots of Γk (x) are exactly those roots of h0 (x) verifying: h1 (x) = 0, . . . , hk−1 (x) = 0, hk (x) = 0 and that if the degree of this greatest common divisor is k then it agrees with Hk (α, y), independently of the value of α. ✷ It is clear from the previous proposition that Φ0 (x) = Γ1 (x) · Γ2 (x) · · · · · Γd−1 (x). This decomposition of the squarefree part of h0 (x) will be called the decomposition with respect to the polynomials h1 (x), h2 (x), . . . , hd−1 (x). Since many of the polynomials Γk (x) are going to be equal to 1, the decomposition of the squarefree part of h0 (x) will be written until the last non constant factor, usually denoted by Γt (x) (1
t < d). Note that since the degree of f with respect to y is d then the maximum posible degree for the (α, y) (with α ∈ C) is d − 1 (those α’s are given by the roots greatest common divisor of f (α, y) and ∂f ∂y of Γd−1 (x)) but since hd−1 (x) is (by definition) dlcofy (f ) then Φd−1 (x) is equal to 1 or to a polynomial without real roots.

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

727

Fig. 3. Generic position.

Next the definition of curve in generic position is presented. It will be the central argument allowing to compute in a very efficient way the topology of the considered real algebraic plane curve. It modifies slightly the definition presented in (Gonzalez-Vega and El Kahoui, 1996) in order to obtain better computational results. Definition 3.4. Let f be a squarefree polynomial in R[x, y]. The real algebraic plane curve defined by f , C(f ), is in generic position if the following two conditions are verified: *1 The leading coefficient of f with respect to y (which is a polynomial in R[x]) has no real roots. *2 For every α ∈ R the number of distinct roots of ∂f (α, y) = 0 f (α, y) = 0, ∂y is 0 or 1. Geometrically, the first condition forbids the presence of vertical asymptotes and the second condition does not allow the existence of two different critical points (vertical tangents, singular points, etc.) sharing the same (real) x-coordinate. In Fig. 3, the first curve is a typical example of non generic position since there are two real roots of the discriminant such that the number of critical points in the vertical lines over these roots is bigger than one. The second curve in Fig. 3 is clearly in generic position since no couple of critical points share the same x-coordinate. Note that second curve is obtained after rotating π/2 radians the first one. The reason why the study of real algebraic plane curves in generic position is much easier than the general case is presented in the next theorem: under such assumption the y-coordinate of a critical point depends rationally on its x-coordinate and if this x-coordinate is a real number then the y-coordinate is real too. Note that this fact is not true if the considered curve is not in generic position. Theorem 3.5. Let f be a squarefree polynomial in R[x, y] such that C(f ) is in generic position. If Hj (x, y) = StHaj (f ) = hj (x)y j + hj,j −1 (x)y j −1 + · · · + hj,0 (x) and (α, β) is a critical point of C(f ) then there exists k such that 1 hk,k−1 (α) . h0 (α) = 0, . . . , hk−1 (α) = 0, hk (α) = 0, β = − k hk (α)

728

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

Proof. It is enough to take into account that if (α, β) is a critical point of C(f ) then α is a root of one of the Γk (x)’s. Moreover this implies (since C(f ) is in generic position) that Hk (α, y) is a polynomial with exactly one root (of multiplicity k). This gives directly the desired description for β. ✷ Corollary 3.6. Let f be a squarefree polynomial in R[x, y] such that C(f ) is in generic position. If α ∈ R verifies that there exists β ∈ C such that f (α, β) = 0,

∂f (α, β) = 0 ∂y

then β ∈ R. Proof. According to Theorem 3.5, there exists k such that β=−

1 hk,k−1 (α) . k hk (α)

Thus, if α ∈ R then β ∈ R.

✷

3.1. Generic position checking The two curves in Fig. 3 can be used to explain why there are only a finite number of bad changes of coordinates with shape: x = x + µy, y = y, such that if C(f ) is not in generic position
then the transformed curve is not in generic position: this  number of bad cases is clearly bounded by 2c where c is the number of distinct critical points. In (Gonzalez-Vega and El Kahoui, 1996) it is shown how algorithmically is always posible to find a change of variable such that the new real algebraic plane curve is in generic position which is briefly reviewed in what follows. Let u be a new variable and g(u, x, y) the polynomial defined by g(u, x, y) = f (x + uy, y). If µ ∈ Z then Cµ (f ) will denote the curve defined by the polynomial g(µ, x, y). Let D(u, x) be the discriminant of g(u, x, y) with respect to the variable y, S(u, x) be its squarefree part and D ◦ (u) the discriminant of S(u, x) with respect to the variable x. In (Gonzalez-Vega and El Kahoui, 1996), the following criteria is presented in order to find µ such that Cµ (f ) is in generic position:

 If µ ∈ Z, degy g(µ, x, y) = deg g(µ, x, y) and D ◦ (µ) = 0 then Cµ (f ) is in generic position. Anyway, since these computations, if not complicated, are quite expensive even for moderate degrees of f (x, y) in terms of computing times a different and more efficient way of checking generic position is proposed here. Next it is shown how running first steps in the algorithm in (Gonzalez-Vega and El Kahoui, 1996) assuming that the given curve is in generic position produces enough information in order to verify if this initial assumption is or not verified.

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

729

Theorem 3.7. Let f (x, y) ∈ R[x, y] be squarefree polynomial of degree d in y such that lcofy (f ) ∈ R (or is a polynomial of positive degree in R[x] without real roots) and Hj (x, y) = StHaj (f ) = hj (x)y j + hj,j −1 (x)y j −1 + · · · + hj,0 (x) (j ∈ {0, . . . , d}) be the Sturm–Habicht sequence of f with respect to y and hd (x), . . . , h0 (x) the corresponding principal Sturm–Habicht coefficients. Let Φ0 (x) be the squarefree part of h0 (x), Φ0 (x) =

h0 (x) , gcd(h0 (x), h0 (x))

and Φ0 (x) = Γ1 (x) · · · · · Γt (x) (1
t < d) be the decomposition of Φ0 (x) with respect to h1 (x), . . . , hd−1 (x). Let α1 < · · · < αs be be the real roots of Γk (x). For the real roots of Φ0 (x) and for every k ∈ {1, . . . , t}, let α1(k) < · · · < αs(k) k k ∈ {1, . . . , t} let rk (x) = −

1 hk,k−1 (x) , k hk (x)

and for k ∈ {1, . . . , t} and j ∈ {1, . . . , sk } let βj(k) = rk (αj(k) ). Then: (1) If f (αj(k) , βj(k) ) = 0 for some k and j then the curve C(f ) is not in generic position. (αj(k), βj(k) ) = 0 for some k and j then the curve C(f ) is not in generic position. (2) If ∂f ∂y (3) The curve C(f ) is in generic position if and only if for every αj(k) the polynomial Hk (αj(k), y) has only one distinct root which is βj(k) . Moreover, in this last case, the critical points of C(f ) are t  


  αj(k) , βj(k) : 1
j
sk .

k=1

Proof. The first two items and the first part of the third are easy since if C(f ) is in generic position then every point (αj(k), βj(k) ) verifies 
f αj(k) , βj(k) = 0,

∂f
(k) (k)  α , βj = 0. ∂y j

For the second part of the third item it is enough to take into account that to say that every Hk (αj(k) , y) has only one distinct root is equivalent to the fact that every system 
∂f
(k)  α ,y = 0 f αj(k) , y = 0, ∂y j has exactly one solution and, thus, to the generic position for C(f ). Example 3.8. The curve defined by the polynomial

✷

730

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

 

f (x, y) := 125y 6 − 300xy 5 + −16 + 315x 2 y 4 + 16x − 184x 3 y 3 
+ 63x 4 − 4x 2 y 2 − 12x 5 y + x 6 having the polynomial
 Φ0 (x) = 3645x 6 − 21168x 4 + 20480x 2 − 1024 x = Γ1 (x)Γ3 (x) as the squarefree part of h0 (x) is in generic position because • for any real root α of Γ1 (x) (in this case there are 6 real roots), defining β=−

2α(2835α 6 − 257424α 4 + 323584α 2 − 16384) h1,0 (α) =− , h1 (α) 70875α 6 + 961680α 4 − 1309952α 2 + 66560

the following equality H1 (α, y) = h1 (α)(y − β) is verified; • for any real root α of Γ3 (x) (in this case α = 0), defining β=−

8α(−32 + 543α 2 ) 1 h3,2 (α) =− = 0, 3 h3 (α) 1024 − 12360α 2 + 1575α 4

the following equality H3 (α, y) = h3 (α)(y − β)3 is verified. The topological graph corresponding to this curve can be found in Fig. 8 corresponding to Pol16 . 3.2. Topology computation for curves in generic position First it is shown how regular points over real roots of the discriminant can be characterized and easily computed for curves in generic position. This is done by presenting an uniform set of univariate polynomials without multiple roots and easy to compute, whose real roots are the y-coordinates of the regular points. Definition 3.9. Let α and β be two new variables. For k ∈ {1, . . . , d − 1} the polynomial Fk (α, β, y) is defined as the quotient of dividing f (α, y) by (y − β)k (using y as main variable). The computation of the coefficients of each Fk (α, β, y) requires no more than the iterated application of Horner/Ruffini’s rule as shown by the next proposition. Proposition 3.10. The coefficients of each polynomial Fk (α, β, y) =

d−k  u=0

G(k−1) (α, β)y d−k−u u

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

731

can be computed in an inductive way by the following relations (u = 0, . . . , d − k, v = 1, . . . , k − 1): u u   (v−1) u−l (0) u (α, β) = G (α, β)β and G (α, β) = β + ad−l (α)β l−1 . G(v) u u l l=0

l=1

Proposition 3.11. With the notations of Theorem 3.7 and if C(f ) is in generic position then: (1) If k ∈ {1, . . . , t} and j ∈ {1, . . . , sk } then  Fk αj(k), βj(k) , y


=

f (αj(k), y) (y − βj(k) )k

being this polynomial division exact. (2) The polynomial Fk (αj(k) , βj(k) , y), of degree d − k, has no multiple roots. Let ν(k,j ) be the number of its real roots (0
ν(k,j )
d − k). (3) For every k ∈ {1, . . . , t} and j ∈ {1, . . . , sk }, if ν(k,j ) = 0 then the curve C(f ) has no regular points (k,j ) ) < . . . < γν(k,j be the real roots of the polynomial on the vertical line x = αj(k) . Otherwise, let γ1 (k,j) Fk (αj(k), βj(k) , y). Then the regular points of C(f ) are: t  


(k,j )

αj(k) , γi

 : 1
j
sk , 1
i
ν(k,j ) .

k=1

Proof. Since C(f ) is in generic position then every polynomial f (αj(k), y) has only one multiple root, βj(k) , of multiplicity k. This implies directly that Fk (αj(k) , βj(k) , y) has no multiple roots and gives the description in the third item about the set of regular points on the vertical lines defined by the real roots of the discriminant. ✷ Finally last proposition shows how to perform the branch counting for curves in generic position. Let α1 < · · · < αs be the real roots of the discriminant. For every i ∈ {1, . . . , s − 1}, let bi denote the number of branches of C(f ) in the interval (αi , αi+1 ). In order to consider all the intervals, including (−∞, α1 ) and (αs , +∞), two values are added to the sequence α1 < · · · < αs : α0 = α1 − 1 and αs+1 = αs + 1 (for example) and define b0 and bs accordingly. Proposition 3.12. Suppose that C(f ) is in generic position and consider two consecutive real roots αi < αi+1 (with 1
i
s − 1) of the discriminant of f with respect to y such that αi = αj(k1 1 ) and αi+1 = αj(k2 2 ) (if k1 = k2 then j2 = j1 + 1). Suppose also that on the vertical line x = αi = αj(k1 1 ) there are νk1 ,j1 regular points and on the vertical line x = αi+1 there are νk2 ,j2 regular points. Then: (1) The number of branches to the left and to the right of any regular point is equal to 1. (2) The number of branches to the right of the critical point (αj(k1 1 ) , βj(k1 1 ) ) is equal to bi − ν(k1 ,j1 ) . (3) The number of branches to the left of the critical point (αj(k2 2 ) , βj(k2 2 ) ) is equal to bi − ν(k2 ,j2 ) . Proof. The first item of the proposition is obvious due to the Implicit Function Theorem. bi is equal to the number of real roots of f (δi , y) (as polynomial in y) where δi is any fixed value between αi and

732

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

Fig. 4. Diagram illustrating the proof of Proposition 3.12.

αi+1 . This number is provided by the sign changes of the sequence h0 (δi ), . . . , hd (δi ) (see Theorem 2.5). Using the first item, we obtain that the number of branches to the right of the critical point (αj(k1 1 ) , βj(k1 1 ) ), is equal to bi minus the number of branches (over (αi , αi+1 )) which begin at a regular point on the vertical line x = αj(k1 1 ) , i.e., bi − ν(k1 ,j1 ) . Analogously, the number of branches to the left of the critical point (αj(k2 2 ) , βj(k2 2 ) ) is equal to bi minus the number of branches (over (αi , αi+1 )) which end at a regular point on the vertical line x = αj(k2 2 ) , i.e., bi − ν(k2 ,j2 ) . ✷ Remark that the condition of having at most one critical point on any vertical line allows this extremely easy half-branch counting. The situation to the left of α1 (respectively to the right of αs ) is particularly simple since all the points of C(f ) over α0 (respectively αs ) are regular. Next example clarifies, in a particular case, the proof of Proposition 3.12. Example 3.13. Suppose that νk1 ,j1 , the number of regular points on the vertical line x = αi = αj(k1 1 ) , is 9, νk2 ,j2 = 7, and the position of the critical and regular points on the vertical lines x = αi and x = αi+1 is as follows (see Fig. 4): (k1 ,j1 )

< γ2

(k2 ,j2 )

< γ2

γ1 γ1

(k1 ,j1 )

< γ3

(k1 ,j1 )

< βj(k1 1 ) < γ4

(k2 ,j2 )

< γ3

(k1 ,j1 )

(k2 ,j2 )

< γ4

(k2 ,j2 )

(k2 ,j2 )

< γ5

(k1 ,j1 )

< γ5

(k1 ,j1 )

< γ7

(k2 ,j2 )

< γ7

< γ6

< βj(k2 2 ) < γ6

(k1 ,j1 ) (k2 ,j2 )

(k1 ,j1 )

< γ8

(k1 ,j1 )

< γ9

.

Suppose also that bi , the number of branches in the interval (αi , αi+1 ), is 12. (k ,j ) As the number of branches to the left and to the right of a regular point is 1, (αi , γ9 1 1 ) will (k ,j ) (k ,j ) (k ,j ) be connected with (αi+1 , γ7 2 2 ) and (αi , γ8 1 1 ) will be connected with (αi+1 , γ6 2 2 ). Analogously (k ,j ) (k ,j ) for every l ∈ {1, 2, 3}, (αi , γl 1 1 ) needs to be connected with (αi+1 , γl 2 2 ). The remaining 4 points

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743 (k1 ,j1 )

γ4

(k1 ,j1 )

, . . . , γ7

733

on the vertical line x = αi , which are above βj(k1 1 ) , will be connected with the critical (k ,j )

(k ,j )

point (αi+1 , βj(k2 2 ) ). Also the remaining 2 points (αi+1 , γ4 2 2 ) and (αi+1 , γ5 2 2 ) on the vertical line x = αi+1 , which are under βj(k2 2 ) , will be connected with the critical point (αi , βj(k1 1 ) ). So far 11 branches in the interval (αi , αi+1 ) have been generated. As bi = 12, the remaining branch will connect the critical points (αi , βj(k1 1 ) ) and (αi+1 , βj(k2 2 ) ). So the number of branches which start at βj(k1 1 ) (i.e., the number of branches to the right of the critical point (αi , βj(k1 1 ) )) is equal to

(k ,j ) (k ,j )  2 corresponding to γ4 2 2 and γ5 2 2 + 1 corresponding to βj(k2 2 ) = 3 = 12 − 9 = bi − νk1 ,j1 . And the number of branches which end at βj(k2 2 ) (i.e., the number of branches to the left of the critical point (αi+1 , βj(k2 2 ) )) is equal to

(k ,j ) (k ,j )  4 corresponding to γ4 1 1 , . . . , γ7 1 1 + 1 corresponding to βj(k1 1 ) = 5 = 12 − 7 = bi − νk2 ,j2 .

4. The algorithm The algorithms computing the topology of the curve defined in an implicit way by a bivariate squarefree polynomial f (x, y) usually require the solving of polynomial equations of type f (α, y) = 0, where α is a real root of the discriminant, and thus with singular solutions. Moreover the determination of what happens to the right and to the left of every critical point requires the computation and manipulation of the real roots of f (α, y) = 0 with α a real root of the discriminant: if using the symbolic approach this is extremely expensive in terms of computing times and if using a purely numerical approach this is a source of stability problems and inacuraccies since one or more real roots of f (α, y) = 0 are singular. Our algorithm avoids the two problems before mentioned by initially performing a linear change of coordinates such that the number of critical points on every vertical line is bounded by 1: i.e., the considered curve is in generic position. If this condition is verified then • every β, for (α, β) a critical point, can be described in an explicit way as a rational function of the corresponding α (see 3.5), • the polynomial f (α, y) is replaced (in a symbolic and explicit way) by a polynomial without multiple roots (see 3.11), and • the behaviour of the curve to the right and to the left of every vertical line x = α is easily determined by knowing what happens in every open interval determined by two consecutive real roots of the discriminant. It is important to quote that the computations producing the critical points of f allow to check if the given curve is in generic position (see Theorem 3.7). In this way the time consuming computations described in (Gonzalez-Vega and El Kahoui, 1996) to produce a change of coordinates such that the transformed curve is indeed in generic position are avoided.

734

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

4.1. The description of the algorithm TOP Next the algorithm TOP computing the topology of C(f ) is described. Input: f ∈ Z[x, y] a squarefree polynomial. •1 Change of variable (making f monic with respect to y if it is not the case). *1 Let µ ∈ Z such that f (x + µy, y) = τy d + ad−1 (x)y d−1 + · · · + a1 (x)y + a0 (x) with τ = 0. Note that there is only a finite number of bad choices for µ. The polynomial f (x + µy, y) will be denoted again by f (x, y). •2 Sturm–Habicht sequence of f and discriminant computation. *1 Let d be the y-degree of f , Hd (x, y), . . . , H0 (x, y) Hj (x, y) = StHaj (f ) = hj (x)y j + hj,j −1 (x)y j −1 + · · · + hj,0 (x) the Sturm–Habicht sequence of f with respect to y (computed by using Subresultant algorithm; see (Lombardi et al., 2000) or (Gonzalez-Vega et al., 1990, 1994)) and hd (x), . . ., h0 (x) the corresponding principal Sturm–Habicht coefficients. *2 Let Φ0 (x) be the squarefree part of h0 (x) (which is the discriminant of f with respect y). •3 Generic position test and critical point description. *1 Let Φ0 (x) = Γ1 (x) · · · · · Γt (x) (1
t
d − 1) be the decomposition of Φ0 (x) with respect h1 (x), . . . , hd−1 (x). be the *2 Let α1 < · · · < αs be the real roots of Φ0 (x). For every k ∈ {1, . . . , t}, let α1(k) < · · · < αs(k) k real roots of Φ0 (x) which are real roots of Γk (x). *3 For k ∈ {1, . . . , t} let rk (x) = −

1 hk,k−1 (x) . k hk (x)

For k ∈ {1, . . . , t} and j ∈ {1, . . . , sk } let βj(k) = rk (αj(k)). *4 If 
∂f
(k) (k)  α , βj = 0 f αj(k) , βj(k) = 0 or ∂y j for some k and j then the curve is not in generic position and the algorithm is restarted (from the step [•1 ]) after applying the change of variable x = x − y, y = y, with f (x − y, y) denoted again by f (x, y), being already monic with respect to y. *5 Otherwise:

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

735

*5.1 If for every k ∈ {1, . . . , t} and j ∈ {1, . . . , sk } the following condition is verified: 


k Hk αj(k) , y = hk αj(k) y − βj(k) then the curve C(f ) is in generic position and the critical points of C(f ) are: t  


  αj(k), βj(k) : 1
j
sk .

k=1

*5.2 Otherwise the curve is not in generic position and the algorithm is restarted (from the step [•1 ]) after applying the change of variable x = x − y, y = y, with f (x − y, y) denoted again by f (x, y), being already monic with respect to y. •4 Regular point computation. *1 For every k ∈ {1, . . . , t} and j ∈ {1, . . . , sk }, let Fk (αj(k), βj(k) , y) be the polynomial obtained by dividing f (αj(k), y) and (y − βj(k) )k . These polynomials have only simple roots. The coefficients of these polynomials d−k  f (αj(k) , y) 

(k) (k)  (k) d−k−u = G(k−1) αj , βj βj Fk αj(k) , βj(k) , y = u (k) k (y − βj ) u=0

can be computed in an inductive way by the following relations (u = 0, . . . , d − k, v = 1, . . . , k − 1): G(v) u




 αj(k), βj(k) =

u 


(k) (k)  (k) u−l G(v−1) αj , βj βj l

and

l=0

G(0) u




 αj(k) , βj(k)

u = βj(k)

+

u 


 l−1 ad−l αj(k) βj(k) .

l=1 (k,j )

*2 For every k ∈ {1, . . . , t} and j ∈ {1, . . . , sk }, let γ1

) < · · · < γν(k,j be the real roots of the (k,j)

polynomial Fk (αj(k), βj(k) , y). The regular points of C(f ) on the vertical lines defined by the real roots of the discriminant are: t  


(k,j ) 

αj(k), γi

 : 1
j
sk , 1
i
ν(k,j ) .

k=1

•5 Branch number computation. *1 Two values are added to the sequence α1 < · · · < αs (in order to consider all the intervals, including (−∞, α1 ) and (αs , +∞)): α0 = α1 − 1 and αs+1 = αs + 1. For every i ∈ {0, . . . , s}, bi represents the number of branches over the interval (αi , αi+1 ). *2 The number of branches to the right of the critical point (αj(k) , βj(k) ) is equal to bi − ν(k,j ) , where αi = αj(k) .

736

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

*3 The number of branches to the left of the critical point (αj(k), βj(k) ) is equal to bi−1 − ν(k,j ) , where αi = αj(k) . (k,j ) *4 The number of branches to the left and to the right of any regular point (αj(k), γi ) is equal to 1. Output: Once the critical points, regular points and number of branches in each interval are known, the graph characterizing the topological structure of the curve C(f ) has been determined: the knots of the graph are the critical points and the regular points computed in [•3 ] and [•4 ]; the process of connecting them by edges is completely determined by the number of branches in each interval through the use of Proposition 3.12 as presented in [•5 ] and illustrated in Example 3.13 and Fig. 4. The algorithm terminates after a finite number of changes of variables since there are only a finite number of bad configurations (i.e., those such that Cµ (f ) is not in generic position: two distinct critical points with the same x-coordinate) for the set of critical points. 4.2. A detailed example Let f (x, y) = y 4 − 6xy 2 + x 2 − 4x 2 y 2 + 24x 3 . Next the Sturm–Habicht sequence is displayed together with the principal Sturm–Habicht coefficients:
 h4 (x) = 1. H4 (x, y) = StHa4 (P ) := y 4 − 6x + 4x 2 y 2 + 24x 3 + x 2 ;
 3 2 h3 (x) = 4. H3 (x, y) = StHa3 (P ) := 4y − 12x + 8x y;
 h2 (x) = 48x + 32x 2 . H2 (x, y) = StHa2 (P ) := 48x + 32x 2 y 2 − 16x 2 − 384x 3 ;
 H1 (x, y) = StHa1 (P ) := 1536x 3 − 1280x 4 − 768x 5 + 512x 6 y; h1 (x) = 1536x 3 − 1280x 4 − 768x 5 + 512x 6 . H0 (x, y) = StHa0 (P ) := 98304x 11 − 585728x 10 + 1253376x 9 − 1126400x 8 + 344064x 7 + 16384x 6 = h0 (x). The discriminant of f with respect to y agrees with stha0 (P ) and its squarefree part is
 Φ0 (x) = 24x 3 − 71x 2 + 45x + 2 x. The decomposition of Φ0 (x) with respect to the principal Sturm–Habicht coefficients produce the polynomials Γ1 (x), Γ2 (x), Γ3 (x) as follows: Γ1 (x) = 24x + 1,

Γ2 (x) = 2 − 3x + x 2 ,

Γ3 (x) = x.

The real roots of Φ0 (x) are α1(1) = −0.0416666,

α1(2) = 1.0,

α2(2) = 2.0,

α1(3) = 0.0

and next the candidates to be the y-coordinates of the critical points are determined by using the formulae in Theorem 3.7: β1(1) = 0.0,

β1(2) = 0.0,

β2(2) = 0.0,

β1(3) = 0.0.

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

737

Now the generic position checking is made: 


 H1 α1(1), y = −0.1148700y = h1 α1(1) y − β1(1) but




2 H2 α1(2), y = 80.0y 2 − 400.0 = h2 α1(2) y − β1(2) = 80.0y 2 .

So C(f ) is not in generic position. The process is restarted after applying the change of variable x = x − y. The new polynomial f (x − y, y) will be denoted again by f (x, y): f (x, y) = −3y 4 + 66xy 2 − 18y 3 + x 2 − 2xy + y 2 − 4x 2 y 2 + 8y 3 x + 24x 3 − 72x 2 y. The Sturm–Habicht sequence together with the principal Sturm–Habicht coefficients are:

 H4 (x, y) = −3y 4 + (8x − 18)y 3 + −4x 2 + 66x + 1 y 2 + −72x 2 − 2x y + 24x 3 + x 2 ; h4 (x) = −3.
 h3 (x) = −12. H3 (x, y) = −12y 3 + (24x − 54)y 2 + −8x 2 + 132x + 2 y − 72x 2 − 2x;
 2
 2 3 2 H2 (x, y) = −288x − 2160x − 2988 y + 192x + 4176x + 7296x + 108 y h2 (x) = −288x 2 − 2160x − 2988. − 1728x 3 − 3984x 2 − 108x;  H1 (x, y) = −1536x 6 − 163584x 5 − 849408x 4 − 1548864x 3 − 997152x 2 − 50112x − 2016 y


+ 96768x 6 + 584832x 5 + 1241856x 4 + 915360x 3 + 44064x 2 + 2016x; h1 (x) = −1536x 6 − 163584x 5 − 849408x 4 − 1548864x 3 − 997152x 2 − 50112x − 2016. H0 (x, y) = −294912x 11 + 1314816x 10 + 4423680x 9 − 7698432x 8 − 28071936x 7 − 19228416x 6 − 801792x 5 − 32256x 4 ; h0 (x) = H0 (x, y). The discriminant of f with respect to y agrees with stha0 (P ) and its squarefree part is
 Φ0 (x) = −x 384x 7 − 1712x 6 − 5760x 5 + 10024x 4 + 36552x 3 + 25037x 2 + 1044x + 42 . The decomposition of Φ0 (x) with respect to the principal Sturm–Habicht coefficients does not produce new factors and thus Φ0 (x) = Γ1 (x). The real roots of Φ0 (x) are α1(1) = −1.318767,

α2(1) = 0.0,

α3(1) = 3.255278,

α4(1) = 5.711823

and next the candidates to be the y-coordinates of the critical points are determined by using the formulae in Theorem 3.7: β1(1) = −2.189490,

β2(1) = 0.0,

β3(1) = 2.262675,

β4(1) = 3.699158.

The verification that C(f ) is in generic position is done in the following way (see Theorem 3.7): 


 H1 α1(1), y = −42482.773331y − 93015.636774 = h1 α1(1) y − β1(1) , 


 H1 α2(1), y = −2016y = h1 α2(1) y − β2(1) , 


 H1 α3(1), y = −221168480.633628y + 500432584.768754 = h1 α3(1) y − β3(1) , 


 H1 α4(1), y = −2273405902.547792y + 8409688762.652116 = h1 α4(1) y − β4(1) . Next the polynomials F1 (αj(1), βj(1) , y) giving the y-coordinates of the regular points over the real roots of the discriminant are computed:

738

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743


 F1 α1(1) , β1(1) , y 
F1 α2(1) , β2(1) , y 
F1 α3(1) , β3(1) , y 
F1 α4(1) , β4(1) , y

= −3y 2 − 15.413194y − 11.119532, = −3y 2 − 18y + 1, = −3y 2 − 5.533827y + 163.777613, = −3y 2 + 5.499633y + 329.219980, (1,j )

together with its real roots γi γ1(1,1) = −4.269618, γ1(1,2) = −6.055050, γ1(1,3) = −8.368322, γ1(1,4) = −9.599104,

(i.e., the y-coordinates of the regular points):

γ2(1,1) = −0.868112, γ2(1,2) = 0.055050, γ2(1,3) = 6.523713, γ2(1,4) = 11.432315.

Thus the searched critical points of C(f ) are   (−1.318767, −2.189490), (0.0, 0.0), (3.255278, 2.262675), (5.711823, 3.699158) and the regular points of C(f ) are  (−1.318767, −4.269618), (−1.318767, −0.868112), (0.0, −6.055050), (0.0, 0.055050),

 (3.255278, −8.368322), (3.255278, 6.523713), (5.711823, −9.599104), (5.711823, 11.432315) .

Finally, by using the polynomials hj (x) and one interior point inside every interval −∞ = α0 < α1(1) < α2(1) < α3(1) < α4(1) < α5 = +∞, it is obtained that the number of branches over every interval is 2, 4, 4, 2 and 4, respectively. For example, taking δ1 =

α1(1) + α2(1) = −0.659383 2

(α1(1) < δ1 < α2(1) ), it is obtained that h3 (δ1 ) = −12, h4 (δ1 ) = −3, h0 (δ1 ) = −321086.868773

h2 (δ1 ) = −1688.949900,

h1 (δ1 ) = −98782.931186,

which implies directly, applying Theorem 2.5 and Definition 2.4, that the number of branches over the interval (α1(1) , α2(1)) is 4 (i.e., 4 sign permanences minus 0 sign variations). All this information produces the graph displayed in Fig. 5, corresponding to Pol4 .

5. Implementation and experimentation The algorithm TOP depends only, first, on the computation of the Sturm–Habicht sequence of f (x, y) with respect to y, second, on the decomposition of Φ0 (x) (squarefree part of the discriminant of f (x, y) with respect to y) with respect to the principal Sturm–Habicht coefficients of f (x, y) and, third, on the computation of the real roots of the polynomials Φ0 (x) and f (α, y) (with α a real root of Φ0 (x)). The algorithm TOP has been implemented in the Computer Algebra System Maple (version 7). It starts

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

739

always with a squarefree polynomial f (x, y) in Q[x, y] (in fact in Z[x, y]) and, due to efficiency reasons, it performs most of the computations by using floating point arithmetic. The only computations which are performed symbolically are: (1) the computation of the Sturm–Habicht sequence of f , and (2) the computation of the squarefree part of the discriminant of f , and its decomposition with respect to the principal Sturm–Habicht coefficients of f . The remaining computations consist in solving numerically different univariate polynomial equations (without multiple roots) or evaluating at these roots some of the polynomials symbolically computed. Initially the chosen precision is 10 digits (the one by default in Maple) and the Maple function fsolve was used to solve the squarefree univariate polynomial equations before mentioned. Once fsolve returns two real roots very close or some numerical evaluation returns some non guaranteed value, the precision is increased by 5 more digits and those computations are performed again (see Example 5.1 below). The experimental results obtained with this implementation are very good. Table 1 contains the results of the experimentation of the Maple implementation on the polynomials in Table 4. The meaning of the different columns in Table 1 is the following: • • • •

Time—means the computing time in seconds, d—denotes the y-degree of the polynomial defining the curve, r—denotes the number of real roots of the discriminant once the curve is in generic position, GP—indicates if the initial polynomial produces a curve in generic position (y: yes) or if the change of variable was required (n: no), and • Precision—denotes the required precision needed to get the desired topological information.

Table 1 Experimental results Pol1 Pol2 Pol3 Pol4 Pol5 Pol6 Pol7 Pol8 Pol9 Pol10 Pol11 Pol12 Pol13 Pol14 Pol15 Pol16

Time

d

r

GP

Precision

2.780 0.110 0.149 0.060 0.039 1.170 0.090 0.569 0.820 0.280 3.190 197.690 0.041 3.780 0.030 0.200

8 5 5 4 4 8 4 6 8 6 6 8 8 6 6 6

28 4 5 4 3 8 2 6 13 5 4 4 3 8 3 7

y y y n n n n y n n y n y n y n

15 10 15 10 10 15 10 20 15 10 15 40 10 25 10 10

740

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743 Table 2 Computing time comparison Pol4 Pol5 Pol6 Pol9 Pol11 Pol15

Time TOP

Time in (Hong, 1996)

0.030 0.039 1.170 0.820 3.190 0.030

0.067 0.117 0.917 0.733 6.267 0.050

Table 3 Two real roots of the discriminant very close α2 and α3

Digits 10 20 30 40

0.9090909091 · 10−8

0.9090909091 · 10−8 0.90909090909090909091 · 10−8 0.90909090909090909091 · 10−8 0.909090909090909090909090909091 · 10−8 0.909090909090909090909090909091 · 10−8 0.9090909090909090909090909090909073334909 · 10−8 0.9090909090909090909090909090909108483272 · 10−8

Fig. 5. Graphs corresponding to Poli , i = 1, 2, 3, 4.

Fig. 6. Graphs corresponding to Poli , i = 5, 6, 7, 8.

Table 2 presents a comparison between the computing times obtained by using the algorithm TOP and by applying the method in (Hong, 1996). It is worth to remark that algorithm TOP presents in general a better behaviour. This is specially relevant if moreover it is taken into account that the computing time of TOP includes the verification that the considered curve is or is not in generic position.

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

741

Fig. 7. Graphs corresponding to Poli , i = 9, 10, 11, 12.

Fig. 8. Graphs corresponding to Poli , i = 13, 14, 15, 16.

The considered examples are presented explicitly in Table 4. Finally Figs. 5, 6, 7 and 8 contain the graphs generated by Maple corresponding to the polynomials in Table 4. Example 5.1. From Table 1, the reader will notice the relatively large computing time corresponding to the polynomial 1 y. 10000000 The reason is found in the fact that in order to separate the real roots of the discriminant it was necessary to go until a precision of 40 digits: note that this polynomial of degree 58 and with exactly four different real roots is known to be squarefree but the function fsolve returns two equal floating point real numbers as real roots until the precision of 40 is attained (see Table 3). Pol12 = y 8 + 25y 4 x 4 + 61y 6 + 62x 5 y 4 + x 6 − 5x 2 y 5 − 11xy +

6. Conclusions The new algorithm presented in this paper has shown a very good practical performance in terms of computing times taking into account that the obtained results are free of numerical errors despite the use of numerical methods to deal with the algebraic numbers involved in the computations. It is worth to remark that it works perfectly on curves with singularities. Two are currently the main weakness of the algorithm in the presented version. The first one is the computation of the real roots of the discriminant: if the degree of f is d then the degree of the discriminant of f with respect to y is bounded by d 2 . Thus for a bivariate polynomial of degree 20, the discriminant is a univariate polynomial whose degree could be up to 400. The computation of the real roots of the discriminant through a matricial formulation as a generalized eigenvalue problem to overcome this weakness is under investigation. The second weakness is the requirement that the initial

742

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

Table 4 Polynomials used in the experimentation Index 1

Polynomial 2 + 7x − 14x 3 + 7x 5 − x 7 − (7 − 42x 2 + 35x 4 − 7x 6 )y − (16 + 42x − 70x 3 + 21x 5 )y 2 − (−14 + 70x 2 − 35x 4 )y 3 − (−20 − 35x + 35x 3 )y 4 − (7 − 21x 2 )y 5 − (8 + 7x)y 6 + y 7 + y 8

2

y 5 + (−x − 1)y 4 + (−2x + 1)y 3 + (2x + 1)y 2 + (2x − 1)y − 2x − 1

3

y 5 + (−x − 1)y 4 + (−2x + 1)y 3 + (2x + 1)y 2 + (2x − 1)y − x − 1

4

y 4 − 6y 2 x + x 2 − 4y 2 x 2 + 24x 3

5

x 4 + 2y 2 x 2 − x 2 + y 4 − 4y 2

6

−3 + 12y 2 + 2y 4 − 12y 6 + y 8 + 12x 2 − 28y 2 x 2 + 12y 4 x 2 + 4y 6 x 2 − 18x 4 + 20y 2 x 4 + 2y 4 x 4 + 12x 6 − 4x 6 y 2 − 3x 8

7

x 5 − y 3 x 2 + xy 4 − x 4 + y 3 x − y 4 + x 3 + yx 2 + y 2 x − 2x 2 − y 2 + x

8

23039 4 359 6 2 4 279756x + 279936xy 4 − 559692y 2 x − 15583x 3 + 217x 5 + 130218 5 x − 10 x + 10 x + 370656y 37333439 2 4 6 4 2 2 3 − 3726432 y 2 − 72774y 2 x 2 + 12947 5 5 y x + 1296y + 46728y x + 15588y x + 100

9

x 8 + 4x 6 y 2 + 6y 4 x 4 + 4y 6 x 2 + y 8 − 4x 6 − 12y 2 x 4 − 12y 4 x 2 − 4y 6 + 16y 2 x 2

10

x 6 + y 2 x 4 − y 4 x 2 − 2x 4 − y 6 + 2y 4 + x 2 − y 2 + xy

11

−2 − 5x − 3y − 6xy 4 − 2y 2 x + y 4 x 4 + 7x 3 y 4 + 6x 3 − x 5 − 6x 2 − 6x 4 − 7x 6 + 4y 3 + 4y 4 + 5y 5 + 3y 2 − 3x 6 y 2 + 3xy + 2y 2 x 2 − 5yx 3 − 2x 6 y − 7x 5 y 3 − 7x 5 y 2 − 4xy 5 − y 3 x 3 + 4yx 2 − 3y 3 x 2 − 6yx 5 − 7y 2 x 4 + y 6 − 3y 4 x 2 − y 3 x + 5yx 4 − 6y 2 x 3 − 3x 6 y 3 + 5x 6 y 4 − 3x 2 y 5 + 4x 5 y 4 − 6x 6 y 5 − 4x 5 y 5 + 2x 4 y 5 + 5x 3 y 5 + 2y 3 x 4

12

1 y 8 + 25y 4 x 4 + 61y 6 + 62x 5 y 4 + x 6 − 5x 2 y 5 − 11xy + 10000000 y

13

y 8 − xy + x 2

14

−2x 6 y 3 + x 5 y 4 + x 7 y + y 4 x 4 + y 6 x 2 − 2y 3 x 4 − 2x 2 y 5 − xy 6 + y 2 x 4 − y 4 x 2 − y 2 x 3 − yx 3 + y 4 − yx 2

15

−x 3 + x 2 y − x 2 − xy 2 + x + y 3 − y 2 − y + 1

16

x 6 + 3y 2 x 4 + 3y 4 x 2 + y 6 − 4y 2 x 2

polynomial f must have integer/rational coefficients while, usually, in problems coming from practice coefficients are floating point real numbers: again the previous matrix formulation can be very useful to overcome this difficulty. The application of this algorithm to two concrete problems in Computer Aided Geometric Design is under consideration: the first one is the computation of the topology of an arrangement of a finite family of planar curves (as in (Keyser et al., 2000)) and the second one is the design of an algorithm allowing a fast characterization of the topology for the sectioning of a NURBS (Non-Uniform Rational B-Spline Surface).

L. Gonzalez-Vega, I. Necula / Computer Aided Geometric Design 19 (2002) 719–743

743

Acknowledgements The authors are very grateful to the two referees for their extraordinary detailed reports, including numerous relevant suggestions, on the original submitted manuscript.

References Arnborg, S., Feng, H., 1988. Algebraic decomposition of regular curves. J. Symbolic Comput. 1–2, 131–140. Arnon, D., McCallum, S., 1988. A polynomial time algorithm for the topological type of a real algebraic curve. J. Symbolic Comput. 5, 213–236. Bajaj, C., Hoffmann, C.M., Lynch, R.E., Hopcroft, J.E.H., 1988. Tracing surface intersections. Computer Aided Geometric Design 5, 285–307. Cellini, P., Gianni, P., Traverso, C., 1991. Algorithms for the shape of semialgebraic sets: a new approach. In: Lecture Notes in Computer Science, Vol. 539. Springer-Verlag, pp. 1–18. Cucker, F., Gonzalez-Vega, L., Rosello, F., 1991. On algorithms for real algebraic plane curves. In: Effective Methods in Algebraic Geometry. In: Progress in Math., Vol. 94. Birkhausser, pp. 63–88. Feng, H., 1992. Decomposition and computation of the topology of plane real algebraic curves, PhD Thesis. The Royal Institute of Technology, Stockholm, Sweden. Gianni, P., Traverso, C., 1983. Shape determination of real curves and surfaces. Ann. Univ. Ferrara Sez VII Sec. Math. XXIX, 87–109. Gonzalez-Vega, L., El Kahoui, M., 1996. An improved upper complexity bound for the topology computation of a real algebraic plane curve. J. Complexity 12, 527–544. Gonzalez-Vega, L., Lombardi, H., Recio, T., Roy, M.-F., 1990. Specialisation de la suite de Sturm et sous-resultants (I). Informatique Theorique et Applications 24, 561–588. Gonzalez-Vega, L., Lombardi, H., Recio, T., Roy, M.-F., 1994. Specialisation de la suite de Sturm et sous-resultants (II). Informatique Theorique et Applications 28, 1–24. Gonzalez-Vega, L., Lombardi, H., Recio, T., Roy, M.-F., 1998. Determinants and real roots of univariate polynomials. In: Caviness, B., Johnson, J. (Eds.), Quantifier Elimination and Cylindrical Algebraic Decomposition. In: Texts and Monographs in Symbolic Computation. Springer-Verlag, pp. 300–316. Hong, H., 1996. An efficient method for analyzing the topology of plane real algebraic curves. Math. Comput. Simulation 42 (4– 6), 571–582. Keyser, J., Culver, T., Manocha, D., Krishnan, S., 2000. Efficient and exact manipulation of algebraic points and curves. Computer-Aided Design 32 (11), 649–662. Keyser, J., Krishnan, S., Manocha, D., 1999a. Efficient and accurate b-rep generation of low degree sculptured solids using exact arithmetic: I—representations. Computer Aided Geometric Design 16 (9), 841–859. Keyser, J., Krishnan, S., Manocha, D., 1999b. Efficient and accurate b-rep generation of low degree sculptured solids using exact arithmetic: II—computation. Computer Aided Geometric Design 16 (9), 861–882. Lombardi, H., Roy, M.-F., Safey el Din, M., 2000. New structure theorem for subresultants. J. Symbolic Comput. 29 (4–5), 663–690. Loos, R., 1982. Generalized polynomial remainder sequences. Computer Algebra (Computing Suplementum) 4, 115–138. Roy, M.-F., 1996. Basic algorithms in real algebraic geometry and their complexity: from Sturm’s theorem to the existential theory of reals. In: Lectures in Real Geometry. In: de Gruyter Exp. Math., Vol. 23. de Gruyter, pp. 1–67. Roy, M.-F., Szpirglas, A., 1990. Complexity of the computation of cylindrical decomposition and topology of real algebraic curves using Thom’s Lemma. In: Lecture Notes in Math., Vol. 1420. Springer-Verlag, pp. 223–236. Sakkalis, T., 1991. The topological configuration of a real algebraic curve. Bull. Australian Math. Soc. 43 (1), 37–50.