- Email: [email protected]
The Cayley method in computer aided geometric design
Computer Aided Geometric Design 1 (1984) 309-326 North-Holland
309
The Cayley method in computer aided geometric design Yves DE MONTAUDOUIN and Wayne TILLER Structural Dynamics Research Corporation, 2000 Eastman Drive, Milfora~ OH 45150, U.S.A. Received 18 December 1984 Revised 9 April 1985 Abstract. A thorough treatment of the Cayley method of elimination is given, with particular emphasis on its implementation to solve problems in computational geometry. An easily programmable formula for the Cayley matrix is derived, it is shown how to extract the birational correspondence from the Cayley matrix, numerous examples are worked, and the connection between the Cayley method and the Euclidean Algorithm is explained.
Keywords. Elimination, Cayley method, birational correspondence, Euclidean Algorithm, polynomials, intersection problems, implicitization.
1. Introduction
We assume the reader is familiar with basic algebraic concepts such as tings, fields, unique factorization domains (UFD), polynomials in one and more variables, and the greatest common divisor (gcd) of two polynomials. We use the following notation: R: an arbitrary UFD (this includes rings of polynomials in one or more variables over an arbitrary field), K: an arbitrary field, C: the field of complex numbers, K[x]: the ring of polynomials in the variable x with coefficients in the field K, K{x}:the field of rational functions in the variable x with coefficients in the field K, K [ x l .... , xn]: the ring of polynomials in the n variables xa,..., xn with coefficients in K. We may also view this ring as the ring of polynomials in the one variable xn with coefficients in the ring K [ x 1 , . . . , X n _ l ] . Hence, we can write K [ x i , . . . . x=] = K [ x 1.... , x,_ll[x,]. If P ~ K[x], deg(P) denotes the degree of P. If P ~ K [ x , y], deg(P) is the total degree of P, degx(P ) is the partial degree in x, and degy(P) the partial degree in y. Algebraic elimination is a classical method of eliminating a common variable from two polynomial equations. The PUrpose of this paper is to examine the Cayley method of elimination, particularly from the viewpoint of its use and implementation in a geometric modeling environment. Most geometric modeling systems use either polynomial or rational functions to represent curves and surfaces. In such systems most geometric problems can be formulated as systems of polynomial equations. Elimination can then be used to either solve the system or to reduce them in some manner. We list some typical examples and, where appropriate, reference previous works which have used elimination to solve these problems. 0167-8396/84/$3.00 © 1984, Elsevier Science Publishers B.V. (North-Holland)
Y. De Montaudouin, W. Tiller / The Cayley method
310
(1) Line~surface intersection Let S(u, v) = (x(u, v), y(u, v), z(u, v)) be a surface represented by x-, y- and z-coordinate functions which are polynomials in the parameters u and v. Represent the line L as the intersection of two planes whose equations are F(x, y, z ) = alx + blY + clz-t- d I = 0 and G(x, y, z) = a2x + b2y + CEZ + d2 = 0. Then if we substitute x(u, v), y(u, v) and z(u, v) into F(x, y, z) and G(x, y, z), we obtain two polynomials in u and v: (S)
V(u,v)=O,
a(u,v)=O.
The solutions to this system are the parameter values (ui, vi) corresponding to the points of intersection of L and S(u, v). Elimination can be used to solve (S). One of the variables, say v, is eliminated from (S), yielding a polynomial R ~ C[u] whose roots are the u-solutions to (S). Subsequently, the set of v-solutions are found and properly paired with the u-solutions to give the solutions to (S). There are two classical methods for determining the v-solutions once the u-solutions have been found: (a) Let u l , . . . , un be the roots of R(u). For each u i, compute the gcd Di(v) of P ( u i , v) and 1 . , vki(i) be the roots of Di(v). Q(ui, v). Since u i is a root of R(u), k(i) = deg(Di) > 0. Let vi,.. Then (u~, Vij'~k(i) J]=l are solutions to (S), and the set (u~, v/), 1 <~i <~n, 1 < j <~k ( i ) is the complete solution set of (S). Kajiya [Kajiya '82] used this method to solve the line/surface intersection problem. H e used the Euclidean Algorithm for computing the gcd of P(u i, v) and Q(ui, v). (b) The other method is to simply repeat the elimination, this time eliminating u to get T(v) ~ C[v]. After finding the roots of T(v), we have a set of u and v values which must be properly paired to yield the solutions to (S). This can be done by substituting all possible pairs (ui, Vl) into (S) and calculating the residues. This method is more stable than the gcd method, but it has the disadvantage that it is difficult to decide " W h a t is a solution?" i.e. " W h e n are the residues zero?" One contribution of this paper is to introduce a third method which we believe is algebraically and computationally superior to the other two methods. The key lies in the following observation: if we are eliminating v, then the v-solutions corresponding to the roots of R(u) can be expressed as ratios of certain sub-determinants of R(u) (R(u) is obtained by expanding a determinant whose elements are polynomials in u). Hence, if we eliminate v we get R(u) and v(u) = f ( u ) / g ( u ) , where f , g ~ C[u]; and if we eliminate u we get T(v) and u(v) = p ( v ) / q ( v ) , with p, q ~ C[v]. These expressions result from the well-known birational correspondence, which we will discuss in more detail below.
(2) Implicitization of a curve Let C(t) = (x(t), y(t)) be a curve whose x- and y-coordinate functions are polynomials in the parameter t. Then we can set F(x, y, t ) = x - x ( t ) and G(x, y, t ) = y - y ( t ) . Then F, G ~ C[x, y, t], and the system F(x, y, t ) = O, G(x, y,, t ) = 0 represents the curve C(t). If we eliminate t from the system we get H(x, y ) = 0, where H E C[x, y]. H is the implicit equation of the curve. The birational correspondence can be extracted during the elimination, and it will have the form t(x, y ) = f ( x , y ) / g ( x , y), with f, g ~ C[x, y]. This problem has been considered in the classical literature, and more recently in the C A D / C A M literature [Sederberg et al. '84, '85].
(3) Curve/curve intersection Let the two xy-planar curves be denoted by Cl(U)=(Xl(U ), Yl(U)) and CE(V)= (xz(v), yE(v)). There are two ways elimination can be used to intersect C 1 and Cz:
:II. De Montaudouin, W. Tiller / The Cayley method,
311
(a) Implicitize one of the curves, say C1. This yields H(x, y) and u(x, y) = f ( x , y)/g(x, y). Then substitute x2(v ) and Y2(o) into H(x, y); the result is a polynomial P(v). Find the roots (vi} of P ( v ) which lie in the valid parameter range of C2. Use these (vi} and x2(v), y2(v) to calculate the points ((x i, Yi)). Substitute these points into the equation u = f / g to obtain the set ( u i ). Let { uj ) be the subset of the { ui ) which lie in the valid parameter range of C 1. Then the ((xj, yj)) corresponding to the ((uj, vj.)) are the desired intersection points. (b) The coordinate functions can be considered as polynomials in two variables, i.e. xl(u, v), yl(u, v), x2(u , v), y2(u, V) E C[u, v], and degv(xl) = dego(y~) = degu(x2) = degu(y2) = 0. Then form the system:
P(u, v ) = x l ( u , v ) - X z ( U , v)=O
and
Q(u, v)=y~(u, v ) - y 2 ( u , v)=O.
This is equivalent to the system (S) above, and it can be solved in the same manner, i.e. eliminate v to get R(u) and v ( u ) = f ( u ) / g ( u ) . In the literature we have only found references to method (a), but we implemented both methods and found (b) to be more efficient.
(4) Surface/surface intersection Let the two surfaces be parametrically defined by Sl(S, t) = (Xl(S, t), yl(s, t), Zl(S , t)) and S2 ( u, v) = ( x : ( u, v ), yz ( U, v ), z 2( u, v )), where the coordinate functions are all polynomials. As in the curve/curve example above, there are two ways to apply elimination to determine the curve of intersection of $1 and $2. Analogous to method (b), we can write:
P(s, t, u, v ) = x l ( s , t ) - x 2 ( u , v ) = o , Q(s, t, u, v ) = y l ( s , t ) - y 2 ( u , v ) = O ,
R(s, t, u, v)=z,(s, t)-z2(u, v)=o. We can then eliminate u from P and Q to get F(s, t, v), and u from Q and R to get G(s, t,'v). Then eliminate v from F and G to get H(s, t ) = 0. This equation contains (as a subset of its real solutions)the points in the parameter space of St whose image is the curve of intersection between $1 and $2. By eliminating s and t from P, Q and R we could obtain the corresponding equation in the parameter space of S 2. See Fig. 1.
Sds,~(u,v) t) Hls,tl=Oj =S
v~'/Klu,v)=O -LI
Fig. 1. The curve of intersection of two parametric surfaces has ;epresentations H(s, t)= 0 and K(u, v ) = 0 in the respective parameter planes of the two surfaces.
312
Y. De Montaudouin, W. Tiller / The Cayley method
(5) Project point to surface Let P = (x, y, z) be a point and S(u, v) a surface whose coordinates are parametric polynomials. Let D(u, v) denote the square of the distance from P to the points of S(u, v). Denote by Du(u, v) and Do(u , v) the partial derivatives of D(u, v) with respect to u and v respectively. Then the points of S(u, v) which are locally a minimum distance from P are among the solutions of the system given by D,,(u, v) = Do(u, v) = 0. This is equivalent to system (S) above.
2. The
Cayley
method
Let R be an arbitrary U F D and P, Q two polynomials in R[x]. The problem in elimination is to set up the so-called resultant of P and Q, denoted by Res(P, Q). Res(P, Q) is an expression in the coefficients of P and Q (i.e. elements of R), and the vanishing of Res(P, Q) is a necessary and sufficient condition for P and Q to have a gcd of degree > 0. The two best k n o w n methods of forming the resultant are Sylvester's and Cayley's methods. In both methods the resultant is calculated as a determinant. Sylvester's method is better suited for theoretical expositions, and a thorough treatment of it can be found in many books [Jacobson '74, Lelong-Ferrand & Amaudies '78, Uspensky '48]. The Cayley method is preferable for a computer implementation because it yields a smaller determinant. However, thorough descriptions of the Cayley method seem to be lacking in both the mathematics and C A D / C A M literature. The best exposition we are aware of is in [Ldong-Ferrand & Amaudies '78]. Our aim here is to make the Cayley method more accessible to the C A D / C A M implementor by presenting a complete, self-contained and as elementary as possible treatment of it. In particular, we: - give a detailed derivation of th6 formula used to calculate the dements of the Cayley matrix; this formula can be easily programmed, - show how to extract the birational correspondence from sub-determinants of the Cayley matrix (and show how to use it to solve geometric problems), - carefully explain special cases such as deg(gcd)> 1, - prove pertinent theorriffs~ - give m a n y simple examples to clarify the concepts, - discuss the connection between the Cayley method and the Euclidean Algorithm for finding the gcd. Before attacking the general problem, we will start with a simpler system. Let P, Q ~ C[x], i.e. polynomials in one variable with complex coefficients. We now show how the Cayley method is used to determine if- P and Q have common roots, and if so, to calculate them. Assume that d e g ( P ) = d e g ( Q ) = m (this assumption simplifies the discussion, it will be explained at the end of this section). Then P ( x ) = Z,k=OakXmk and Q(x)-- Z-,l=0~'~m blxt. Let x 0 ~ C. Then f ( x ) = P ( x o ) Q ( x ) - e ( x ) Q ( x o ) is a polynomial of degree at most m, and f ( x 0 ) = 0. Hence, g(x) = f ( x ) / x - Xo) is a polynomial in x of degree at most m - 1, and the symmetry of f i x ) implies that the coefficients of g(x) are (m - 1)-th degree polynomials in x 0. So we can write:
g(x)=
m~l (m-1 . ) j=o
E cj,x ,
i=0
/
where the cj~ are functions of the a k and bI. N o w assume that x o is a common root of P ( x )
Y. De Montaudouin, W. Tiller / The Cayley method
313
and Q(x). Then g ( x ) is zero for all x, and it follows that m--1
cjix~ = O for all j , O <~j <~m - 1 . i=O
This yields the homogeneous system of linear equations:
c00
C0mlI0 ]
•
"
tc, _l,0
...
•
~ 0 .
Cm_l,m_lJk'X 21
The system has a nontrivial solution because of the 1 in the unknowns column, ar;d it follows that det[cji] = 0, where det denotes the determinant. [cji] is called the Cayley matrix. The following theorem gives an even stronger result than the one just derived.
Theorem 1. Let P, Q ~ C[x] and assume that deg(P) = deg(Q) = m > 0. Let [cji ] be the Cayley matrix. Then the following two conditions are equivalent: (1) P and Q have r common roots (distinct or not). (2) The rank of [cji ] is m - r. Proof. We skip the proof here, as we will prove the theorem below for the more general case: P,Q ~ R[x], for R an arbitrary U F D . Example 1. Let P ( x ) = - 1 + x 2 and Q ( x ) = 2 - 3x + x 2. Then
f(x)= P(xo)Q(x) - P(x)a(xo) = ( - 1 + x o 2 ) ( 1 - 3x + x 2 ) - ( - 1 = ( 3 x 0 - 3x2) + ( - 3
+ x 2 ) ( 2 - 3x o + x 2)
+ 3X2o)X + ( 3 -
3x0)x 2.
Dividing by x - x o gives g ( x ) = ( - 3 + 3x0) + (3 - 3x0)x. Since g ( x ) = 0 for all x, it follows that (-= 3 + 3x0) = (3 - 3 x 0 ) = 0, or equivalently:
A=[
-3
3 / is the Cayley matrix, and clearly r a n k ( A ) = 1. From 3 - 3 x 0 = 0 , 3 - 3 .I that x 0 = 1 is the common root of P ( x ) and Q(x).
it follows
1.
We will now d.erive the Cayley matrix for the general case: P, Q ~ R[x], where R is an arbitrary U F D . There is one subtle nuisance here, namely the concept of a root. Although it is not at all obvious for general R what a root of P ~ R[x] is, we will continue to use the term. This is mathematically justifiable, and in fact, the use of such symbolic roots is a common and powerful practice in algebra. There are approaches which do not rely on the concept of a root, but they are less elegant, less intuitive, and mathematically more sophisticated. Again, we assume that d e g ( P ) - d e g ( Q ) = m > 0 (see below)• Then P ( x ) = Ek=oakX m k and Q ( x ) =~.t=obtx, m l where a k, b t ~ R. As above, we fix Xo and form the functions f ( x ) = P ( x o ) Q ( x ) - P ( x ) Q ( x o ) and g ( x ) = f ( x ) / ( x - Xo). We will get the Cayley matrix from g(x). F r o m the expressions for P and Q, we obtain: m,m
P(xo)Q(x)=
~. k,l=O,O
m,m
akbtxkx t and
P(x)Q(xo)=
E k,l=O,O
atbk x k x '
Y. De Montaudouin, IrK.Tiller / The Cayley method
314
and hence m~D'l
f(x)=
(akbt-atbk)X~X'.
E k,l=0,0
Regrouping coefficients, we can sum on the set 0 ~< k < l ~< m, and obtain: f(x)=
(akbt--albk)(xgxt--x~ xk)
E O~
E (akbl--atbk)xkxk(xl-k--xl-k) O<~k
•
A n d using the algebraic identity: n-1
x n - Y " = ( x - - y ) ( x n - I + x " - 2 Y + "'" +Y"-I)=(x--Y) E xn-l--iY i, i=0
yields
g(x)=f(x)/(X-Xo) l--k-1
=
Z
(akbt-atbk)Xg xk ~-, xt-k-l-ix~
O~k
i=0
£ E O<~k<~m-1 k+l<~l
l-k--1 E (akbt--atbk) xk+ixt-l-i. i=0
N o w the goal is to obtain a polynomial in x, so we need to reverse the index i so that the powers of x are increasing, and we need to move the /-summation to the outside. Setting I = l - 1 - i yields: 1--1
g(x) =
E "E £ (akbl--albk)xlo -1-I+kxI" O~
A n d noting the equivalence of the summation sets:
O<~k~m-l} (!~I<.m-1 k + l <~l <~m ** <~k <~I k <~I <~l - 1 +l~
g(x)=
I
m
E
Z
Z
I=0
k=0
l=I+1
(akbt--atbk)Xto -1-1+kxz"
N o w if x 0 is a c o m m o n root of P and Q, then zero, and hence: I
m
E
E
(akbt-atbg) xt-l-I+k=O,
g(x) is a polynomial in x which is identically forallI, O ~ l ~ m - 1 .
k=O 1=I+ 1
With the goal of obtaining a polynomial in x0, we set L = 1 - 1 - 1 + k and get: 1 m--1 - l + k
Y'~
Z
k=O
L=k
(akbr+l+,-k-- a r+l+I-kbk)XOL =0"
Noting the equivalence of the summation sets:
O~k~I k<~L<~m-l-I+k
} ** { O , L , m - 1 max(0, L - m + l + I ) < k < ~ r f f m ( L ,
I)
315
Y. De Montaudouirg W. Tiller/ The Cayley method
a n d setting m n -- m i n ( L , I ) and m x = max(0, L - m + 1 + I ) , it follows that: m--1
inn
E E (akbL+a+I--k L=0 k=mx
-
-
a z + l + z - k b k ) x oi. = O.
N o w set inn
ClL =
E
(akbL+l+Z-k--az.+l+I-kbk).
k= mx
T h e n [ciL ] is the Cayley matrix. N o t e that [czr] is symmetric, i.e. cn. = czj. If x 0 is a c o m m o n r o o t of P a n d Q, then we have the m × m system of h o m o g e n e o u s equations: m--1
Y" czLx~=O,
O~I~m-1.
L=O
W e n o w state the f u n d a m e n t a l theorem for the Cayley matrix: T h e o r e m 2. Let P, Q ~ R [ x ] , R an arbitrary UFD. Assume that d e g ( P ) = d e g ( Q ) --- m > 0. Let [czL ] be the Cayley matrix. Then the following two conditions are equivalent: (1) The degree of the gcd of P and Q is r. (2) The rank of [cir] is m - r. Before proving T h e o r e m 2, we give an example.. E x a m p l e 2. T a k e P, Q ~ C[x][y]:
3 P(x, y)=(x2-x)+(-5x+7)y+(3x-6)y
2+y3=
y, akyk=O,
k=0 3
a(x, y)= (-2x
+ 2)+(3x-
1)y +(-7x
+ 10)y 2 +y3=
E bty t = O ,
1=0 where a o = x 2 -
b o = - 2 x + 2,
x,
a 1 = - 5x + 7,
b 1 = 3x - 1,
a 2 = 3x - 6,
b 2 = - 7x + 10,
a3=1,
b 3 = 1.
N o w the 3 × 3 C a y l e y matrix is given by: a0b I - alb o
aob2-a2bo
a°b3 - d3b° ]
aob2 - aEb o (aob3-a3bo)+(albE-aEbl)
aob 3 - a3b o
axb 3 - a3b 1 J ' a2b 3 - a3b 2
aab3 - a3b 1
a n d hence,
M(x) =
3x 3 - 14x 2 + 25x - 14 - 7 x 3 + 23x 2 - 28x + 12
xEwx--2
- T x 3 + 23x 2 - 28x + 12 27x 2 -8x+
_
77x + 62
8
E x p a n d i n g d e t ( M ( x ) ) , we get the polynomial: F ( x ) -- 16 456 - 64 2 8 4 x + 107 794x 2 - - 102 085x 3 q _ 59 995x 4
- 2 2 285x 5 + 4899x 6 - 490x 7.
x2+x--2] -8x+ 8 1. l O x - 16 J
Y. De Montaudouin, W. Tiller / The Cayley method
316
F(x) not identically zero implies the matrix M(x) is nonsingular, and from Theorem 2 it follows that the degree of the gcd of P and Q is zero, i.e. P and Q have no common components. N o w x = 1 is a root of F, and M(1)
!2 • 0 Since r a n k ( M ( 1 ) ) = 2, it follows that P(1, y) and Q(1, y), as polynomials in C[y], have one common root. Indeed, =
P(1, y ) = 2 y - 3y 2 + y 3 = y ( y _ 2 ) ( y - 1), Q(1, y ) = 2 y + 3 7 2 + y 3 =y(y + 2 ) ( y + 1), and the common root is y = 0. Hence, (1, 0) is a solution to the system P, 0. Also, x = 2 is a root of F, and M(2) =
[4_, i]
-8 16 • 4 -8 Since rank(M(2))= 1, it follows that P(2, y) and Q(2, y) have two common roots. Since P(2, y ) = 2 - 3y + y3 __ ( y _ 1)2(y + 2), 0 ( 2 ; y ) = - 2 + 5 y - 4y 2 + y 3 = ( y _ 1 ) 2 ( y _ 2), it follows that y = 1 is a double root. Hence, (2, 1) is a double solution. By computing the tangent lines to P and Q at (2, 1), it can be seen that they touch tangent to one another at this point. Fig. 2 shows the two curves P(x, y) and Q(x, y). P has two real branches; Q has one. For Q one can determine that
3- y +o in the vicinity of infinity. This shows that the line x = y/7 + 73/49 is an asymptote to Q, and it also gives the position of Q relative to the asymptote. The two branches of P are asymptotic to the line
y=
+2- y
and the parabola 16 11( x = - 3y 2 + - ~ y - ~ -
52) 243y "
Q(x,y)
l'>J ,
~
=,.x
Fig. 2. Two intersecting algebraic curves:
Plx;y)
P(x, y ) = ( x 2 - x ) + ( - 5x + 7 ) y 4-(3x - 6 ) y 2 4- y3,
Q(x, y ) = ( - 2 x 4- 2)+(3x - 1 ) y + ( - 7 x + 1 0 ) y 2 + y3.
Y. De Montaudouin, W. Tiller / The Cayley method
317
We now turn our attention to the proof of Theorem 2. We are not aware of any elementary proofs of this theorem in the literature. A recent proof can be found in [Griffiths '81], but because of the different context, his proof is more complicated than ours. Our proof requires a n u m b e r of propositions and lemmas; propositions will be stated without proof, either because they are well-known from algebra, or because they are easy to prove; lemmas will be proved.
Definition. A polynomial P ( x I , . : . , Xn) is called symmetric in its n variables if any permutation of the variables leaves the polynomial unchanged.
Definition. A polynomial P ( x l , . . . , xn) is called homogeneous of degree m if P ( ) ~ x l , . . . , X x , ) = ~mp(x I .... , Xn).
For the rest of this section we use the following notation: n
P, Q~R[x],
P(x)=
E ai xi,
n
Q(x)=
i=0
Y~ bjx j, j=O
and a , , b, not equal to zero. We denote by a~,..., a n the roots of P, and fla . . . . , fin the roots of Q. We will interchangeably denote the Cayley matrix of P and Q by [ClL], C ( P , Q), or C ( a o , . . . , an, bo . . . . , bn).
Proposition 3. P ( x ) = an(x - ax) - - • ( x - an), and n
a n - 1 = -- E
an an-2 = an
Oti,
i=1 E
otiotj,
l <~i
n
a---9-°= ( -- 1)n 1-I ai = ( -- 1)no~1o~2 • • • ~n" an i=1
Proposition 4. The expressions for a k / a , , 0 <~k <~n - 1, are symmetric, homogeneous polynomials in the roots a i. Furthermore, any symmetric, homogeneous polynomial in the roots a i has exactly one representation as a polynomial in the n indeterminants ( a n _ l / a n , . . . , ao/an). Conversely, assume that F ( a 0 , . . . , an) is a polynomial in the ai which is homogeneous of degree p. Then F ( a o , . . . , a n ) = a ~ f ( a o / a n , . . . , a n - 1 / a , ) . And we can replace the a j a n with their expressions in the at, and we obtain: aP,f ( ao/an, . . . , a , _ l / a , ) = aP,g( a l , . . . , an). We will apply these techniques to the Cayley matrix.
Lemma 5. (1) det(C(P, Q)) is homogeneous of degree n separately in a o , . . . , a n and bo,... , b,. (2) I f the expressions in the a i and flj are substituted for the a i and bj in d e t ( C ( P , Q)), the result is a homogeneous polynomial of degree n 2. In shorthand, det(C(ao,...,an,
bo,...,bn))
= a~b~ d e t ( C ' ( a o / a n , . . . , a , _ l / a n , b o / b n , . . . , bn_t/bn) ) = a'~b~' det( C"( al,. .., a , , fla . . . . . fin)).
Y. De Montaudouin, W. Tiller / The Cayley method
318
Proof. (1) follows easily from the definition of the Cayley matrix:
XP(x)~O.(Xo) -- Xo
= X # ( e ( x ) e ( x o ) -- Xo e(xo)Q(X
)
implies C(XP, # Q ) = hl~C(P , Q), and it follows that det(C(TkP, i~Q))=det(X#C(P, Q ) ) = xn~n d e t ( C ( P , Q)). N o w the proof of (2). For any k, 0 ~
ak/a . will be of degree n - k. The same holds for fli and bk/b n. Let QL be an dement of C(P, Q). Then by definition, n--1
CIL = E
akbl+L+l--k-- aI+L+l-kbk.
k=0
If we replace the a~, bi (0 <~ i, j ~< n - 1) with their roots a~,/3i (and the factors a , , bn), then we obtain for ClL a polynomial of degree (n - k ) + (n - (1 + L + 1 - k)) = 2n - 1 - ( 1 + L) in the a i and flj. The expansion of d e t ( C ( P , Q)) has the form: (*)
det(C(P, a))=
E
sg(°)Co, o(o)Cl,oo)"'" c,-1,o(,-1),
where S,, is the group of all permutations of the integers 0 , . . . , n - 1. N o w choose an arbitrary summand of (*), i.e. fix o ~ S,. Furthermore, we drop sg(o), since it is always equal to + 1 or 1. To complete the proof we need only show that -
deg(c0,,,(0)cl,o(1) ... c,,_ 1 , o ( n - 1)) = n2Since the degree of cz~ is 2n - 1 - ( I + L ) , it follows that deg(C,,o(,)) = 2 n - 1 - 1 - o ( I ) , and hence n--1
deg(c0,o(0)---Cn_l.o(n_l))
=
~ (2n-- 1--1--o(1)) I=0
*
n--1
=n(2n-1)-
n--1
E I-
Z o(1).
I=0
I=0
But Y'.z=oO(1) n-1 n--1 does not depend on o and is equal to Ei=0I. Hence, n--1
deg(Co,o
c._l,o(._l))=n(2n
-
1)-2
I I=0
----2n2-- n-- 2( n ( n -=- n1 2) ")2
[]
Lemma 6. C(P, Q) is singular if and only if P and Q have a common root. Proof. Assume x 0 is a common root. Then from the derivation of C(P, Q), we have the n × n system of homogeneous equations: n--1
~., CILX~=O,
O<.I<~n--1.
L=O
Since x ° = 1, we have a nontrivial solution which implies C(P, Q) is singular. Now assume
319
Y. De Montaudouirg tE. Tiller / The Cayley method
C(P, Q) is singular, i.e. det(C(P, Q ) ) = 0. Since a~ = flj implies det(C(P, Q ) ) = 0, it follows that ai - flj divides det(C(P, Q)). Furthermore, the ( a i - flj), 1 <~i, j <~n, are relatively prime in R [ a l , . . . , an, ill,-.., fin], and hence l-li~j=l(ai- flj) divides det(C(P, Q)). Since deg(I-li,j(ai flj)) = deg(det(c(P, Q)) = n 2, it follows from Lemma 5 that -
0 = det(C(P, Q))=pa~b~X--[(ai-flj), i,j where p is a nonzero integer constant. Therefore, a i = flj for some i, j.
[]
Proposition 7. I f P1,..., Pn are n polynomials which form a system of rank r in the basis 1, x .... , x ,, and if D is any nonzero polynomial, then DP1,.. ., DP n also forms a system of rank r. We now prove the fundamental theorem, Theorem 2: the rank loss of the Cayley matrix is equal to the degree of the gcd. Proof. Let P, Q ~ R[x], d e g ( P ) = d e g ( Q ) = n. Let D ~ R[x] be the gcd of P and Q, and deg(D) = r >I 0. Then P = DP 1, Q = DQ1, and deg(P 0 = deg(Q 0 = n - r. Using the above notation (derivation of the Cayley matrix), we see that D(xo)Pl(Xo)D(x)al(x ) - D(x)Pl(x)D(xo)al(xo) gp,e(x,
Xo) =
x - ~o
= D(x)D(xo)gP1,QI(X , Xo) n--r--1 n--r--1
=D(x)D(xo)
E
E
I=0 L=O
CILX~XI'
where [cIL ] is the Cayley matrix for P1, Q1- Lemma 6 implies that rank([CIL])=n--r. Multiplying through by D ( x ) and D(xo), we obtain
r(rx g,,,Q(~, Xo)-- ,~oE LE0=
)
ClLD(xo)~k D(x)~'.
So [C•L ] is also the Cayley matrix for P and Q in the basis D ( x ) , D ( x ) x .... , D ( x ) x n - ' - l . Rewrite gp,Q(X, XO) as gp,Q(X, XO)= [ D ( x ) , ' " , D ( X ) x n - r - 1 ] [ C I L ]
FID(x0). I
L D(xo)x~) - ' - 1 From Proposition 7 it follows that [ D ( x ) , . . . , D ( x ) x "-'-1] is independent and can be extended to a basis of the polynomials of degree n - 1:
~,Q(x, Xo)--1
fx(x0) = [1,/l(x),...,/,_l(x), D(x), D(x)x,...,D(x)x'-r-q[ei,]
fr-l(XO) D(Xo)
D(xo)xU "-1
Y. De Montaudouin,W. Tiller/ The Cayleymethod
320
where [elL] is the extension of [cIL] obtained by adding r null rows and columns. Clearly, rank([Si~ ]) = rank([cIL]) = n -- r. By a change of basis, we can find an n × n lower triangular invertible matrix A such that 1
D(x)
=
n--1
D(x)xn-,-] and it follows that
gp,Q(X, Xo)=[1, x, .... Xn-1]AT[cqL]AIlXo[ .
Hence, AT[a~L]A is the Cayley matrix for P and Q, and the observation that rank(AT[a~L]A) = rank([a~,]) = n - r completes the proof. []
The following corollary will be of fundamental importance in extracting the birational correspondence from the Cayley matrix (next section). We do not know if this corollary exists in the literature. Coronary 8. I f rank(C(P, Q ) ) = n - r , corner of C ( P , Q) is nonsingular.
then the ( n - r ) - t h
order submatrix in the lower right
Proof. I n t h e proof of Theorem 2 above we showed that C ( P , Q ) = A T C ( P 1, Q1)A, where A is lower triangular and invertible, and C(P1, Q1) is the extension of C(P1, Q1) which contains C(P1, Q1) in its lower right corner and zeros otherwise. Since P1 and Qt are relatively prime, it follows that C ( P 1, Qx) is nonsingular, and hence r a n k ( C ( P 1, Q 1 ) ) = rank(C(P1, Q1))= n - r. The following observation from linear algebra will complete the proof: let B be an arbitrary n × n matrix, and denote by Bp the (n - p ) - t h order submatrix in the lower right corner of B. Let C = ATBA, where A is lower triangular and invertible. Then for every p, 1 ~ p ~< n - 1, rank(Cp) = rank(Bp). But this can be seen by noting that, for i , j > p ,
Cij =
i.e. Cp =
E akibklalj = E akibklalj; l <~k,l<.n p+ l <~k,l<~n (aT)pBpap. []
Throughout this section we have assumed d e g ( P ) = deg(Q). We now show how to remove this restriction. Let P, Q ~ R[x], P ( x ) = •x"m " J, with m > n. Let a l , . . . , a m -,i=o a ix i, Q ( x ) = ~'~j=objx be the roots of P, and ill,---, fin the roots of Q. There are two ways to increase deg(Q) from n to m: (a) Regard Q as a polynomial of degree m, with the m - n leading coefficients equal to zero. (b) Multiply a by x " - n . Definition. The resultant of P and Q, denoted by Res(P, Q), is given by • ,.t \
re+n--1
rct n
m
re
amb il=il j----1 H
Y. De Montaudouin. W. Tiller / The Cayley method
321
Note that this definition does not require deg(P)= deg(Q). We will drop ( - 1 ) m+'-l, since it is not relevant to the discussion. It is easy to verify that
n (**)
ra
Res(P, Q ) = b ~ I - [ P ( f l / ) = a ~ , I " I a ( a i ) . j=l
i=l
When eliminating x, Res(P, Q) is actually what we are looking for when we set up the Cayley matrix; and if deg(P)= deg(Q), it is what we get, i.e. det(C(P, Q))= Res(P, Q) (see the proof of Lemma 6). If the degrees are not equal, we use either method (a) or (b) above. These methods yield a matrix C(P, Q) such that det(C(P, Q ) = F . Res(P, Q), F ~ R; that is, we get an extraneous factor. We now derive this factor. Definition. The reciprocal of a polynomial P, denoted by Rec(P), is the polynomial resulting
from reversing the order of the coefficients. If P ( x ) = a o + a m+
...
. . . +am xm, then R e c ( P ) =
+ a o Xm.
Proposition 9. I f ai are the roots of P, then 1 / a i are the roots of Rec(P). Lemma 10. Res(P, Q ) = Res(Rec(P), Rec(Q)). Proof. m
n
Res(Rec( P ), Rec(Q)) = agb~' Y1 I--I ( 1 / a i - 2/flj )
i=1 j=l
= agbg I-m[ f i (
i=1 j~l
flJ - ai ) Oli~j m
n
= ( a ~ b r ) ( a ~ n b r / a ; b r ) I-I I-I ( ¢ j - - a,) i=l
= Res(P, Q).
j=l
[]
Lennna 11. Let P, Q, T ~ R[x]. Then Res(P, QT) = Res(P, Q) Res(P, T). Proof. Follows from (* *).
[]
Now suppose we use method (a). Let Q denote the polynomial obtained by considering Q as a polynomial of degree m. Theorem 12. det(C(P, Q))= a,~-" Res(P, Q).
Proof. Clearly, R e c ( Q ) = x ~ - " Rec(Q). Then det(C(P, Q)) = Res(P, Q) = Res(Rec(P), Rec(Q)) = Res(Rec( P ), x m-" Rec( a )) = Res(Rec( P ), x ~ - " ) Res(Rec( P ), Rec( Q )) (from Lemma 11). From (* *) it follows that Res(Rec(P), x m - ' ) = (Rec(P(O))) m-" = a m-", and from Lemma 10, Res(Rec(P), R e c ( Q ) ) - R e s ( P , Q). Therefore, det(C(P, Q ) ) = amm - - n Res(P, Q). [] Now method (b). Let Q -- x m-'Q.
Y. De Montaudouin, W. Tiller / The Cayley method
322
T h e o r e m 13. d e t ( C ( P , Q ) ) = a on l - - r l Res(P, Q). Proof.
d e t ( C ( P , Q ) ) = R e s ( P , Q) = R e s ( P ,
xm-nQ)
= Res(P, x m-") Res(P, Q)= = a~'-" R e s ( P , Q).
(P(O))m-"Res(P, Q)
[]
So the extraneous factor, F, is just a~ ' - n if we use m e t h o d (b), a n d a,~ - " if we use m e t h o d (a). In a computer implementation, it should be checked which of the ai, i = 0, m, is simpler. T h e n d e t ( C ( P , Q)) must be divided by a 2 - " Again, we turn to Example 2 (see above). Consider P and Q as polynomials in x; i.e.
P, Q ~ C[yl[x]: 2
P(y, x ) = ( 7 y - 6 y 2 - y 3 ) + ( - 1
- 5y+ 3y2)x + x 2 =
E
ak x~,
k=0 1
a(y,
x)=(2-y+
10y2+y3)+(-2+
3 y - 7 y 2 ) x = ~_,btx'. l=0
m = d e g ( P ) = 2, n = deg(Q) = 1. We eliminate x. T h e 2 × 2 Cayley matrix has the form:
aobl - albo aob2- a2bo] aob2 - a2bo alb2- a2b l j" M e t h o d (a):
b 0 = 2 - y + 10y 2 + y 3 , bl = - 2 + 3y - 7y 2, b2 = 0.
D e n o t i n g the Cayley matrix by M ( y ) , we get:
M(Y)= [2_-25yy+_ 32y21 0 y 2 _ yX5y3 3 +
-2 +y - 10y 2 _y3]
2- 3y+7y 2
]
and d e t ( M ( y ) ) = - 1 2 y + 52y 2 - I45y 3 + 211y 4 - 205y 5 + 169y 6 - 70y 7. m--n = a~ = 1, it follows that d e t ( M ( y ) ) - R e s ( P , Q). Since a m M e t h o d (b):
bo = 0 b 1-- 2 - y + b2= -2+
10y 2 + y 3 , 3 y - 7 y 2.
D e n o t i n g the Cayley matrix by M ( y ) , we get: _~t(y) ---- [
14y -- 19y 2 + 78y 3 -- 54y 4 + 4y 5 -I-y 6
[ -14y
+ 33y 2 - 69y 3 + 45y 4 - 7y 5
-- 14y + 33y 2 -- 69y ~ + 45y 4 - 7y 5 ] 8y - 24y 2 + 43y 3 -- 21y 4
J
Y. De Montaudouin, W. Tiller / The Cayley method
323
d e t ( M ( y ) ) = _ 84y 2 + 436y 3 - 1339y 4 + 2399y 5 - 2846y 6 + 2624y 7 - 1709y 8 + 589y 9 - 70y 10 = (7y - 6y 2 + y 3 ) - d e t ( M ( y ) ) = a 0 • Res(P, Q).
3. T h e birational c o r r e s p o n d e n c e The birational correspondence is often stated a n d / o r motivated using rational curves in the xy-plane. In this context it can be stated as follows: let C(t)= (x(t), y(t)) be a parametric curve, where x(t), y ( t ) ~ K ( t ) , K a field. Assume that for each point on the curve there corresponds only o n e value of t. T h e n there exists f, g ~ K [x, y] such that t f(x, y)/g(x, y). In particular, if K is the field of rational numbers, then the coordinates (X(to), y(to)) of a point on the curve are-rational if and only if to is rational. In the next section we will show that this is a simple consequence of the Euclidean Algorithm; in this section we show that it is a result of eliminating t in the implicitization process (see implicitization example in the Introduction), a n d in particular, we derive from the Cayley matfix a formula which gives t as a function of x a n d y. The unit circle centered at the origin is the standard example for the birational correspondence. It is given by: (1)
x(t)=(1-t2)/(1
+ t 2)
y(t)=2t/(1 +t2).
and
Using substitution it is easy to verify that t = y / ( 1 + x). Let us derive this equation from the Cayley matrix. First, we rewrite (1) as:
F(x, y, t ) = ( x -
1)+(x+
1)t2=0,
G(x, y, t) = y - 2t + yt 2 = O. T h e n F, G E K[x, y][t]. Eliminating t, the Cayley matrix is: M(x, y ) = [ ( x - 1 ) ( - 2 ) - ( O ) ( y )
(x-1)(y)-(x+l)(y)]
(x-
1)(y) - (x + 1)(y)
-2y
2x + 2
__
]
( 0 ) ( y ) - (x + 1 ) ( - 2 )
E x p a n d i n g d e t ( M ( x , y)) and setting it equal to zero yields the circle. A n d f r o m -2y
2x + 2
x 2 + y2"= 1, the implicit-equation of
= 0,
it follows that ( - 2 y ) + (2x + 2)t --- 0, and hence, t = 2y/(2x + 2) = y / ( x + 1). N o t e that we could also have derived t = (1 - x ) / y (from the first line of M(x, y)). But at nonsingularities, these are equivalent since
y/(1 + x ) = y 2 / ( y + xy)= ( 1 - x 2 ) / ( y ( 1 = ((1 - x ) ( 1 +
x))/(y(1
+ x))=
+ x)) (1 -
x)/y.
T h e significance of the birational correspondence for the implicitization process is obvious. Suppose we have a parametric curve C(t) and a set of points {(x i, yi)}, 1 ~< i ~ n, and assume we want to k n o w if they lie on C(t), and if so, what are their corresponding t-values. This is a
Y. De Montaudouin, W. Tiller / The Cayley method
324
time consuming process when done directly from the parametric representation. Using Cayley we get the implicit form P(x, y ) = 0 of C(t) and the corresponding t = f ( x , y)/g(x, y). We substitute the (x;, y~) into P(x, y) to check if they lie on C(t), and if so, t = f / g gives their corresponding parameter values. N o w we look at the more general case, where the Cayley matrix may be larger than 2 x 2. For simplicity, we choose P, Q E K[x][y], K the field of rationals. The generalization to K[x a. . . . , x,][y] is obvious. Suppose d e g y ( P ) = degy(Q) = m > 0. Then
Lcm_l,0(x)
---
C..--1,m--,(X)
..--a
N o w suppose x0 is a value such that det[c~j(Xo)]= O. Then by Theorem 1, P(xo, y) and Q(x o, y), viewed as polynomials in K[y], have a gcd of degree = r > 0, where r is the rank loss of [cij(xo)]. Suppose r = 1. By Corollary 8, we let M 1 denote the ( m - 1)-th order nonsingular submatrix in the lower right corner of [cij(Xo)]. Then
Fy0] 1 c ,o(xo)
(2)
~
"
o
~
~
Lyg,-
.
LCm-l,o(Xo)
and we can solve for Y0 using Cramer's rule: det(A) Yo = d e t ( M 1 ) ' where A is the matrix formed by replacing the first column of M 1 with the right hand side of (2). Clearly, Y0 is rational if x0 is rational. Obviously, we cannot use the above method to solve for Yo if r > 1. In this case there is no birational correspondence. For curve implicitization r = deg(gcd) > 1 implies a singular point on the curve; for curve/curve intersection it implies either that the curves touch tangent to one another or that there are several y-solutions corresponding to one x-solution. We now return to Example 2 of the previous section. W h e n we substitute x o = 1 into the Cayley matrix we obtain: M(1) =
[i 0 0] 12 0
0 • -6
Then r a n k ( M ( 1 ) ) = 2. The last two columns are independent, and if we discard the first row ( k = 1), we get: M I = [12 0
O] -6 "
[1~
Y0 =
Hence, 0
6][Yogi [:] _
and Y°=-
I°0
-6
0
° l=0
-6
Y. De Montaudouin, W. Tiller / The Cayley method
325
This shows that (1, 0) is a solution. For x 0 = 2, rank(M(2))= 1, which implies there is no birational correspondence in this case. This is clear since the gcd of P(2, y) and Q(2, y) is quadratic, i.e. it could have complex roots.
4. The connection between the Cayley method and the Euclidean Algorithm
Let t'1, P2 ~ C[x, y]. We may consider/'1 and P2 as belonging to C(x}[y]; i.e. they are polynomials in the variable y with coefficients which are rational functions in x. N o w C( x } is a commutative field, which implies that C(x }[y] is a Euclidean domain. It follows that we have a Euclidean algorithm in C( x }[y]. Furthermore, C( x }[y] a Euclidean domain implies that it is a principal domain, and hence, we have the following result: if P1 and P2 have no common components (factors) in C[x, y], then they are relatively prime in C(x }[y] and there exist U'1, U~ ~ C{x}[y] such that: U~P1+ U~P2= 1. Now if R is a common multiple of the denominators of the coefficients of U{ and U~, then we can multiply through by R and obtain: U1P1 + U2P2= R, where R ~ C[x] and U1, U2, P1, P2 ~ C[x, y]. Now to find R we can do a Euclidean algorithm in C(x}[y]. Assume degy(P2)>/degy(P1) > 0. Dividing P2 by P1 we get: P2 = P1Q1 + R1, Q1, R1 ~ C(x}[y] and degy(R1) < degy(P1). If we repeat this enough times, we eventually get a remainder equal to zero. For illustration, suppose this requires four steps. Then: P2 = P I Q 1 + R1,
P1 = RaQ2 + R2, R 1 = R 2 Q 3 + R 3, R 2 = R 3 Q 4 + O.
Now R 3 is a constant in C{x)[y], that is R 3 ~ C ( x ) , and we have: R 3 = R 1 - R2Q 3
=Ra-Q3(P1-R1Q2 ) = R1 -
Q3PI+
Q302(P2 - P1Qa)
= R 1 -- Q 3 P 1 + Q 3 Q 2 P 2 - P I Q I Q 2 Q 3 = P2 - P I Q 1 - Q 3 P 1 + Q 3 Q 2 P 2 - P 1 Q I Q 2 Q 3
= P I ( - - Q 1 - Q3 - Q 1 Q 2 Q 3 ) + / ' 2 ( 1 + a 2 0 3 ) ,
where Qa, Q2, Q3 ~ c{x}[y] and R 3 ~ C ( x ). Now multiplying through by a common multiple of the denominators of R3, 01, Q2 and Q3 yields: R=P1U1 +P2U2, where U1,
U2 ~ C[x, y] and R ~ C[x]. When we eliminate y from Pa(x, y) and P2(x, y), by the Cayley or any other method, it is R that we find. But R2 is also of interest, as it gives us an expression of y as a function of x. Since d e g y ( R 3 ) = 0, it follows that dege(Rz)>~ 1. Suppose d e g y ( R 2 ) = 1. Then R2(x, y)= Ao(x)+yAa(x), where Ao, A I ~ C ( x } . N o w if Xo is a root of R, and hence, of R3, then R2(xo, y) is the gcd of Pl(xo, y) and P2(x o, y). Hence, if Yo is a root of R2(xo, Y), then (Xo, Yo) is a solution to the system Pl(X, y)=P2(x, y ) = 0 , and Rz(xo, yo)=Ao(xo)+ yAa(xo)--- 0, which implies that y = --Ao(xo)/Aa(xo). This expression is the birational correspondence, and it is equivalent to the expression of y as a ratio of two sub-determinants derived above as a by-product of the Cayley matrix, in the sense that they yield the same value of y for those values of x which are x-solutions to the system Pl(x, y)= P2(x, y ) = 0. The equation y = - A J A x corresponds to a rank loss of one in the Cayley matrix. If there are several
326
Y. De Montaudouin, W. Tiller / The Cayley method
y-solutions for some x-solution, then the gcd R 2 ( x , y) can be of degree k > 1 in y: R z ( x , y ) = A o ( x ) + .. • + y k A k ( x ). This corresponds to a rank loss of k in the Cayley matrix.
Acknowledgement We thank Ron Goldman for reading the original manuscript and making several valuable suggestions, one of which led to the current, stronger version of Corollary 8.
References Griffiths, H.B. (1981), Cayley's version of the resultant of two polynomials, Amer. Math. Monthly 88 (5), 328-338. Jacobson, N. (1974), Basic Algebra 1, Freeman, San Francisco. Kajiya, J.T. (1982), Ray tracing parametric patches, Computer Graphics 16 (3), 245-253. Lelong-Ferrand, J. and Arnaudies, J.M. (1978), Cours de Mathematiques, Tome 1 Algbbre, Dunod, Paris. Sederberg, T., Anderson, D. and Goldman, R. (1984), Implicit representation of parametric curves and surfaces, Computer Vision, Graphics and Image Processing 28, 72-84. Sederberg, T., Anderson, D. and Goldman, R. (1985), Implicitization, inversion, and intersection of rational cubic curves, to appear. Uspensky, J.V. (1948), Theory of Equations, McGraw-Hill, New York.
309
The Cayley method in computer aided geometric design Yves DE MONTAUDOUIN and Wayne TILLER Structural Dynamics Research Corporation, 2000 Eastman Drive, Milfora~ OH 45150, U.S.A. Received 18 December 1984 Revised 9 April 1985 Abstract. A thorough treatment of the Cayley method of elimination is given, with particular emphasis on its implementation to solve problems in computational geometry. An easily programmable formula for the Cayley matrix is derived, it is shown how to extract the birational correspondence from the Cayley matrix, numerous examples are worked, and the connection between the Cayley method and the Euclidean Algorithm is explained.
Keywords. Elimination, Cayley method, birational correspondence, Euclidean Algorithm, polynomials, intersection problems, implicitization.
1. Introduction
We assume the reader is familiar with basic algebraic concepts such as tings, fields, unique factorization domains (UFD), polynomials in one and more variables, and the greatest common divisor (gcd) of two polynomials. We use the following notation: R: an arbitrary UFD (this includes rings of polynomials in one or more variables over an arbitrary field), K: an arbitrary field, C: the field of complex numbers, K[x]: the ring of polynomials in the variable x with coefficients in the field K, K{x}:the field of rational functions in the variable x with coefficients in the field K, K [ x l .... , xn]: the ring of polynomials in the n variables xa,..., xn with coefficients in K. We may also view this ring as the ring of polynomials in the one variable xn with coefficients in the ring K [ x 1 , . . . , X n _ l ] . Hence, we can write K [ x i , . . . . x=] = K [ x 1.... , x,_ll[x,]. If P ~ K[x], deg(P) denotes the degree of P. If P ~ K [ x , y], deg(P) is the total degree of P, degx(P ) is the partial degree in x, and degy(P) the partial degree in y. Algebraic elimination is a classical method of eliminating a common variable from two polynomial equations. The PUrpose of this paper is to examine the Cayley method of elimination, particularly from the viewpoint of its use and implementation in a geometric modeling environment. Most geometric modeling systems use either polynomial or rational functions to represent curves and surfaces. In such systems most geometric problems can be formulated as systems of polynomial equations. Elimination can then be used to either solve the system or to reduce them in some manner. We list some typical examples and, where appropriate, reference previous works which have used elimination to solve these problems. 0167-8396/84/$3.00 © 1984, Elsevier Science Publishers B.V. (North-Holland)
Y. De Montaudouin, W. Tiller / The Cayley method
310
(1) Line~surface intersection Let S(u, v) = (x(u, v), y(u, v), z(u, v)) be a surface represented by x-, y- and z-coordinate functions which are polynomials in the parameters u and v. Represent the line L as the intersection of two planes whose equations are F(x, y, z ) = alx + blY + clz-t- d I = 0 and G(x, y, z) = a2x + b2y + CEZ + d2 = 0. Then if we substitute x(u, v), y(u, v) and z(u, v) into F(x, y, z) and G(x, y, z), we obtain two polynomials in u and v: (S)
V(u,v)=O,
a(u,v)=O.
The solutions to this system are the parameter values (ui, vi) corresponding to the points of intersection of L and S(u, v). Elimination can be used to solve (S). One of the variables, say v, is eliminated from (S), yielding a polynomial R ~ C[u] whose roots are the u-solutions to (S). Subsequently, the set of v-solutions are found and properly paired with the u-solutions to give the solutions to (S). There are two classical methods for determining the v-solutions once the u-solutions have been found: (a) Let u l , . . . , un be the roots of R(u). For each u i, compute the gcd Di(v) of P ( u i , v) and 1 . , vki(i) be the roots of Di(v). Q(ui, v). Since u i is a root of R(u), k(i) = deg(Di) > 0. Let vi,.. Then (u~, Vij'~k(i) J]=l are solutions to (S), and the set (u~, v/), 1 <~i <~n, 1 < j <~k ( i ) is the complete solution set of (S). Kajiya [Kajiya '82] used this method to solve the line/surface intersection problem. H e used the Euclidean Algorithm for computing the gcd of P(u i, v) and Q(ui, v). (b) The other method is to simply repeat the elimination, this time eliminating u to get T(v) ~ C[v]. After finding the roots of T(v), we have a set of u and v values which must be properly paired to yield the solutions to (S). This can be done by substituting all possible pairs (ui, Vl) into (S) and calculating the residues. This method is more stable than the gcd method, but it has the disadvantage that it is difficult to decide " W h a t is a solution?" i.e. " W h e n are the residues zero?" One contribution of this paper is to introduce a third method which we believe is algebraically and computationally superior to the other two methods. The key lies in the following observation: if we are eliminating v, then the v-solutions corresponding to the roots of R(u) can be expressed as ratios of certain sub-determinants of R(u) (R(u) is obtained by expanding a determinant whose elements are polynomials in u). Hence, if we eliminate v we get R(u) and v(u) = f ( u ) / g ( u ) , where f , g ~ C[u]; and if we eliminate u we get T(v) and u(v) = p ( v ) / q ( v ) , with p, q ~ C[v]. These expressions result from the well-known birational correspondence, which we will discuss in more detail below.
(2) Implicitization of a curve Let C(t) = (x(t), y(t)) be a curve whose x- and y-coordinate functions are polynomials in the parameter t. Then we can set F(x, y, t ) = x - x ( t ) and G(x, y, t ) = y - y ( t ) . Then F, G ~ C[x, y, t], and the system F(x, y, t ) = O, G(x, y,, t ) = 0 represents the curve C(t). If we eliminate t from the system we get H(x, y ) = 0, where H E C[x, y]. H is the implicit equation of the curve. The birational correspondence can be extracted during the elimination, and it will have the form t(x, y ) = f ( x , y ) / g ( x , y), with f, g ~ C[x, y]. This problem has been considered in the classical literature, and more recently in the C A D / C A M literature [Sederberg et al. '84, '85].
(3) Curve/curve intersection Let the two xy-planar curves be denoted by Cl(U)=(Xl(U ), Yl(U)) and CE(V)= (xz(v), yE(v)). There are two ways elimination can be used to intersect C 1 and Cz:
:II. De Montaudouin, W. Tiller / The Cayley method,
311
(a) Implicitize one of the curves, say C1. This yields H(x, y) and u(x, y) = f ( x , y)/g(x, y). Then substitute x2(v ) and Y2(o) into H(x, y); the result is a polynomial P(v). Find the roots (vi} of P ( v ) which lie in the valid parameter range of C2. Use these (vi} and x2(v), y2(v) to calculate the points ((x i, Yi)). Substitute these points into the equation u = f / g to obtain the set ( u i ). Let { uj ) be the subset of the { ui ) which lie in the valid parameter range of C 1. Then the ((xj, yj)) corresponding to the ((uj, vj.)) are the desired intersection points. (b) The coordinate functions can be considered as polynomials in two variables, i.e. xl(u, v), yl(u, v), x2(u , v), y2(u, V) E C[u, v], and degv(xl) = dego(y~) = degu(x2) = degu(y2) = 0. Then form the system:
P(u, v ) = x l ( u , v ) - X z ( U , v)=O
and
Q(u, v)=y~(u, v ) - y 2 ( u , v)=O.
This is equivalent to the system (S) above, and it can be solved in the same manner, i.e. eliminate v to get R(u) and v ( u ) = f ( u ) / g ( u ) . In the literature we have only found references to method (a), but we implemented both methods and found (b) to be more efficient.
(4) Surface/surface intersection Let the two surfaces be parametrically defined by Sl(S, t) = (Xl(S, t), yl(s, t), Zl(S , t)) and S2 ( u, v) = ( x : ( u, v ), yz ( U, v ), z 2( u, v )), where the coordinate functions are all polynomials. As in the curve/curve example above, there are two ways to apply elimination to determine the curve of intersection of $1 and $2. Analogous to method (b), we can write:
P(s, t, u, v ) = x l ( s , t ) - x 2 ( u , v ) = o , Q(s, t, u, v ) = y l ( s , t ) - y 2 ( u , v ) = O ,
R(s, t, u, v)=z,(s, t)-z2(u, v)=o. We can then eliminate u from P and Q to get F(s, t, v), and u from Q and R to get G(s, t,'v). Then eliminate v from F and G to get H(s, t ) = 0. This equation contains (as a subset of its real solutions)the points in the parameter space of St whose image is the curve of intersection between $1 and $2. By eliminating s and t from P, Q and R we could obtain the corresponding equation in the parameter space of S 2. See Fig. 1.
Sds,~(u,v) t) Hls,tl=Oj =S
v~'/Klu,v)=O -LI
Fig. 1. The curve of intersection of two parametric surfaces has ;epresentations H(s, t)= 0 and K(u, v ) = 0 in the respective parameter planes of the two surfaces.
312
Y. De Montaudouin, W. Tiller / The Cayley method
(5) Project point to surface Let P = (x, y, z) be a point and S(u, v) a surface whose coordinates are parametric polynomials. Let D(u, v) denote the square of the distance from P to the points of S(u, v). Denote by Du(u, v) and Do(u , v) the partial derivatives of D(u, v) with respect to u and v respectively. Then the points of S(u, v) which are locally a minimum distance from P are among the solutions of the system given by D,,(u, v) = Do(u, v) = 0. This is equivalent to system (S) above.
2. The
Cayley
method
Let R be an arbitrary U F D and P, Q two polynomials in R[x]. The problem in elimination is to set up the so-called resultant of P and Q, denoted by Res(P, Q). Res(P, Q) is an expression in the coefficients of P and Q (i.e. elements of R), and the vanishing of Res(P, Q) is a necessary and sufficient condition for P and Q to have a gcd of degree > 0. The two best k n o w n methods of forming the resultant are Sylvester's and Cayley's methods. In both methods the resultant is calculated as a determinant. Sylvester's method is better suited for theoretical expositions, and a thorough treatment of it can be found in many books [Jacobson '74, Lelong-Ferrand & Amaudies '78, Uspensky '48]. The Cayley method is preferable for a computer implementation because it yields a smaller determinant. However, thorough descriptions of the Cayley method seem to be lacking in both the mathematics and C A D / C A M literature. The best exposition we are aware of is in [Ldong-Ferrand & Amaudies '78]. Our aim here is to make the Cayley method more accessible to the C A D / C A M implementor by presenting a complete, self-contained and as elementary as possible treatment of it. In particular, we: - give a detailed derivation of th6 formula used to calculate the dements of the Cayley matrix; this formula can be easily programmed, - show how to extract the birational correspondence from sub-determinants of the Cayley matrix (and show how to use it to solve geometric problems), - carefully explain special cases such as deg(gcd)> 1, - prove pertinent theorriffs~ - give m a n y simple examples to clarify the concepts, - discuss the connection between the Cayley method and the Euclidean Algorithm for finding the gcd. Before attacking the general problem, we will start with a simpler system. Let P, Q ~ C[x], i.e. polynomials in one variable with complex coefficients. We now show how the Cayley method is used to determine if- P and Q have common roots, and if so, to calculate them. Assume that d e g ( P ) = d e g ( Q ) = m (this assumption simplifies the discussion, it will be explained at the end of this section). Then P ( x ) = Z,k=OakXmk and Q(x)-- Z-,l=0~'~m blxt. Let x 0 ~ C. Then f ( x ) = P ( x o ) Q ( x ) - e ( x ) Q ( x o ) is a polynomial of degree at most m, and f ( x 0 ) = 0. Hence, g(x) = f ( x ) / x - Xo) is a polynomial in x of degree at most m - 1, and the symmetry of f i x ) implies that the coefficients of g(x) are (m - 1)-th degree polynomials in x 0. So we can write:
g(x)=
m~l (m-1 . ) j=o
E cj,x ,
i=0
/
where the cj~ are functions of the a k and bI. N o w assume that x o is a common root of P ( x )
Y. De Montaudouin, W. Tiller / The Cayley method
313
and Q(x). Then g ( x ) is zero for all x, and it follows that m--1
cjix~ = O for all j , O <~j <~m - 1 . i=O
This yields the homogeneous system of linear equations:
c00
C0mlI0 ]
•
"
tc, _l,0
...
•
~ 0 .
Cm_l,m_lJk'X 21
The system has a nontrivial solution because of the 1 in the unknowns column, ar;d it follows that det[cji] = 0, where det denotes the determinant. [cji] is called the Cayley matrix. The following theorem gives an even stronger result than the one just derived.
Theorem 1. Let P, Q ~ C[x] and assume that deg(P) = deg(Q) = m > 0. Let [cji ] be the Cayley matrix. Then the following two conditions are equivalent: (1) P and Q have r common roots (distinct or not). (2) The rank of [cji ] is m - r. Proof. We skip the proof here, as we will prove the theorem below for the more general case: P,Q ~ R[x], for R an arbitrary U F D . Example 1. Let P ( x ) = - 1 + x 2 and Q ( x ) = 2 - 3x + x 2. Then
f(x)= P(xo)Q(x) - P(x)a(xo) = ( - 1 + x o 2 ) ( 1 - 3x + x 2 ) - ( - 1 = ( 3 x 0 - 3x2) + ( - 3
+ x 2 ) ( 2 - 3x o + x 2)
+ 3X2o)X + ( 3 -
3x0)x 2.
Dividing by x - x o gives g ( x ) = ( - 3 + 3x0) + (3 - 3x0)x. Since g ( x ) = 0 for all x, it follows that (-= 3 + 3x0) = (3 - 3 x 0 ) = 0, or equivalently:
A=[
-3
3 / is the Cayley matrix, and clearly r a n k ( A ) = 1. From 3 - 3 x 0 = 0 , 3 - 3 .I that x 0 = 1 is the common root of P ( x ) and Q(x).
it follows
1.
We will now d.erive the Cayley matrix for the general case: P, Q ~ R[x], where R is an arbitrary U F D . There is one subtle nuisance here, namely the concept of a root. Although it is not at all obvious for general R what a root of P ~ R[x] is, we will continue to use the term. This is mathematically justifiable, and in fact, the use of such symbolic roots is a common and powerful practice in algebra. There are approaches which do not rely on the concept of a root, but they are less elegant, less intuitive, and mathematically more sophisticated. Again, we assume that d e g ( P ) - d e g ( Q ) = m > 0 (see below)• Then P ( x ) = Ek=oakX m k and Q ( x ) =~.t=obtx, m l where a k, b t ~ R. As above, we fix Xo and form the functions f ( x ) = P ( x o ) Q ( x ) - P ( x ) Q ( x o ) and g ( x ) = f ( x ) / ( x - Xo). We will get the Cayley matrix from g(x). F r o m the expressions for P and Q, we obtain: m,m
P(xo)Q(x)=
~. k,l=O,O
m,m
akbtxkx t and
P(x)Q(xo)=
E k,l=O,O
atbk x k x '
Y. De Montaudouin, IrK.Tiller / The Cayley method
314
and hence m~D'l
f(x)=
(akbt-atbk)X~X'.
E k,l=0,0
Regrouping coefficients, we can sum on the set 0 ~< k < l ~< m, and obtain: f(x)=
(akbt--albk)(xgxt--x~ xk)
E O~
E (akbl--atbk)xkxk(xl-k--xl-k) O<~k
•
A n d using the algebraic identity: n-1
x n - Y " = ( x - - y ) ( x n - I + x " - 2 Y + "'" +Y"-I)=(x--Y) E xn-l--iY i, i=0
yields
g(x)=f(x)/(X-Xo) l--k-1
=
Z
(akbt-atbk)Xg xk ~-, xt-k-l-ix~
O~k
i=0
£ E O<~k<~m-1 k+l<~l
l-k--1 E (akbt--atbk) xk+ixt-l-i. i=0
N o w the goal is to obtain a polynomial in x, so we need to reverse the index i so that the powers of x are increasing, and we need to move the /-summation to the outside. Setting I = l - 1 - i yields: 1--1
g(x) =
E "E £ (akbl--albk)xlo -1-I+kxI" O~
A n d noting the equivalence of the summation sets:
O<~k~m-l} (!~I<.m-1 k + l <~l <~m ** <~k <~I k <~I <~l - 1 +l~
g(x)=
I
m
E
Z
Z
I=0
k=0
l=I+1
(akbt--atbk)Xto -1-1+kxz"
N o w if x 0 is a c o m m o n root of P and Q, then zero, and hence: I
m
E
E
(akbt-atbg) xt-l-I+k=O,
g(x) is a polynomial in x which is identically forallI, O ~ l ~ m - 1 .
k=O 1=I+ 1
With the goal of obtaining a polynomial in x0, we set L = 1 - 1 - 1 + k and get: 1 m--1 - l + k
Y'~
Z
k=O
L=k
(akbr+l+,-k-- a r+l+I-kbk)XOL =0"
Noting the equivalence of the summation sets:
O~k~I k<~L<~m-l-I+k
} ** { O , L , m - 1 max(0, L - m + l + I ) < k < ~ r f f m ( L ,
I)
315
Y. De Montaudouirg W. Tiller/ The Cayley method
a n d setting m n -- m i n ( L , I ) and m x = max(0, L - m + 1 + I ) , it follows that: m--1
inn
E E (akbL+a+I--k L=0 k=mx
-
-
a z + l + z - k b k ) x oi. = O.
N o w set inn
ClL =
E
(akbL+l+Z-k--az.+l+I-kbk).
k= mx
T h e n [ciL ] is the Cayley matrix. N o t e that [czr] is symmetric, i.e. cn. = czj. If x 0 is a c o m m o n r o o t of P a n d Q, then we have the m × m system of h o m o g e n e o u s equations: m--1
Y" czLx~=O,
O~I~m-1.
L=O
W e n o w state the f u n d a m e n t a l theorem for the Cayley matrix: T h e o r e m 2. Let P, Q ~ R [ x ] , R an arbitrary UFD. Assume that d e g ( P ) = d e g ( Q ) --- m > 0. Let [czL ] be the Cayley matrix. Then the following two conditions are equivalent: (1) The degree of the gcd of P and Q is r. (2) The rank of [cir] is m - r. Before proving T h e o r e m 2, we give an example.. E x a m p l e 2. T a k e P, Q ~ C[x][y]:
3 P(x, y)=(x2-x)+(-5x+7)y+(3x-6)y
2+y3=
y, akyk=O,
k=0 3
a(x, y)= (-2x
+ 2)+(3x-
1)y +(-7x
+ 10)y 2 +y3=
E bty t = O ,
1=0 where a o = x 2 -
b o = - 2 x + 2,
x,
a 1 = - 5x + 7,
b 1 = 3x - 1,
a 2 = 3x - 6,
b 2 = - 7x + 10,
a3=1,
b 3 = 1.
N o w the 3 × 3 C a y l e y matrix is given by: a0b I - alb o
aob2-a2bo
a°b3 - d3b° ]
aob2 - aEb o (aob3-a3bo)+(albE-aEbl)
aob 3 - a3b o
axb 3 - a3b 1 J ' a2b 3 - a3b 2
aab3 - a3b 1
a n d hence,
M(x) =
3x 3 - 14x 2 + 25x - 14 - 7 x 3 + 23x 2 - 28x + 12
xEwx--2
- T x 3 + 23x 2 - 28x + 12 27x 2 -8x+
_
77x + 62
8
E x p a n d i n g d e t ( M ( x ) ) , we get the polynomial: F ( x ) -- 16 456 - 64 2 8 4 x + 107 794x 2 - - 102 085x 3 q _ 59 995x 4
- 2 2 285x 5 + 4899x 6 - 490x 7.
x2+x--2] -8x+ 8 1. l O x - 16 J
Y. De Montaudouin, W. Tiller / The Cayley method
316
F(x) not identically zero implies the matrix M(x) is nonsingular, and from Theorem 2 it follows that the degree of the gcd of P and Q is zero, i.e. P and Q have no common components. N o w x = 1 is a root of F, and M(1)
!2 • 0 Since r a n k ( M ( 1 ) ) = 2, it follows that P(1, y) and Q(1, y), as polynomials in C[y], have one common root. Indeed, =
P(1, y ) = 2 y - 3y 2 + y 3 = y ( y _ 2 ) ( y - 1), Q(1, y ) = 2 y + 3 7 2 + y 3 =y(y + 2 ) ( y + 1), and the common root is y = 0. Hence, (1, 0) is a solution to the system P, 0. Also, x = 2 is a root of F, and M(2) =
[4_, i]
-8 16 • 4 -8 Since rank(M(2))= 1, it follows that P(2, y) and Q(2, y) have two common roots. Since P(2, y ) = 2 - 3y + y3 __ ( y _ 1)2(y + 2), 0 ( 2 ; y ) = - 2 + 5 y - 4y 2 + y 3 = ( y _ 1 ) 2 ( y _ 2), it follows that y = 1 is a double root. Hence, (2, 1) is a double solution. By computing the tangent lines to P and Q at (2, 1), it can be seen that they touch tangent to one another at this point. Fig. 2 shows the two curves P(x, y) and Q(x, y). P has two real branches; Q has one. For Q one can determine that
3- y +o in the vicinity of infinity. This shows that the line x = y/7 + 73/49 is an asymptote to Q, and it also gives the position of Q relative to the asymptote. The two branches of P are asymptotic to the line
y=
+2- y
and the parabola 16 11( x = - 3y 2 + - ~ y - ~ -
52) 243y "
Q(x,y)
l'>J ,
~
=,.x
Fig. 2. Two intersecting algebraic curves:
Plx;y)
P(x, y ) = ( x 2 - x ) + ( - 5x + 7 ) y 4-(3x - 6 ) y 2 4- y3,
Q(x, y ) = ( - 2 x 4- 2)+(3x - 1 ) y + ( - 7 x + 1 0 ) y 2 + y3.
Y. De Montaudouin, W. Tiller / The Cayley method
317
We now turn our attention to the proof of Theorem 2. We are not aware of any elementary proofs of this theorem in the literature. A recent proof can be found in [Griffiths '81], but because of the different context, his proof is more complicated than ours. Our proof requires a n u m b e r of propositions and lemmas; propositions will be stated without proof, either because they are well-known from algebra, or because they are easy to prove; lemmas will be proved.
Definition. A polynomial P ( x I , . : . , Xn) is called symmetric in its n variables if any permutation of the variables leaves the polynomial unchanged.
Definition. A polynomial P ( x l , . . . , xn) is called homogeneous of degree m if P ( ) ~ x l , . . . , X x , ) = ~mp(x I .... , Xn).
For the rest of this section we use the following notation: n
P, Q~R[x],
P(x)=
E ai xi,
n
Q(x)=
i=0
Y~ bjx j, j=O
and a , , b, not equal to zero. We denote by a~,..., a n the roots of P, and fla . . . . , fin the roots of Q. We will interchangeably denote the Cayley matrix of P and Q by [ClL], C ( P , Q), or C ( a o , . . . , an, bo . . . . , bn).
Proposition 3. P ( x ) = an(x - ax) - - • ( x - an), and n
a n - 1 = -- E
an an-2 = an
Oti,
i=1 E
otiotj,
l <~i
n
a---9-°= ( -- 1)n 1-I ai = ( -- 1)no~1o~2 • • • ~n" an i=1
Proposition 4. The expressions for a k / a , , 0 <~k <~n - 1, are symmetric, homogeneous polynomials in the roots a i. Furthermore, any symmetric, homogeneous polynomial in the roots a i has exactly one representation as a polynomial in the n indeterminants ( a n _ l / a n , . . . , ao/an). Conversely, assume that F ( a 0 , . . . , an) is a polynomial in the ai which is homogeneous of degree p. Then F ( a o , . . . , a n ) = a ~ f ( a o / a n , . . . , a n - 1 / a , ) . And we can replace the a j a n with their expressions in the at, and we obtain: aP,f ( ao/an, . . . , a , _ l / a , ) = aP,g( a l , . . . , an). We will apply these techniques to the Cayley matrix.
Lemma 5. (1) det(C(P, Q)) is homogeneous of degree n separately in a o , . . . , a n and bo,... , b,. (2) I f the expressions in the a i and flj are substituted for the a i and bj in d e t ( C ( P , Q)), the result is a homogeneous polynomial of degree n 2. In shorthand, det(C(ao,...,an,
bo,...,bn))
= a~b~ d e t ( C ' ( a o / a n , . . . , a , _ l / a n , b o / b n , . . . , bn_t/bn) ) = a'~b~' det( C"( al,. .., a , , fla . . . . . fin)).
Y. De Montaudouin, W. Tiller / The Cayley method
318
Proof. (1) follows easily from the definition of the Cayley matrix:
XP(x)~O.(Xo) -- Xo
= X # ( e ( x ) e ( x o ) -- Xo e(xo)Q(X
)
implies C(XP, # Q ) = hl~C(P , Q), and it follows that det(C(TkP, i~Q))=det(X#C(P, Q ) ) = xn~n d e t ( C ( P , Q)). N o w the proof of (2). For any k, 0 ~
ak/a . will be of degree n - k. The same holds for fli and bk/b n. Let QL be an dement of C(P, Q). Then by definition, n--1
CIL = E
akbl+L+l--k-- aI+L+l-kbk.
k=0
If we replace the a~, bi (0 <~ i, j ~< n - 1) with their roots a~,/3i (and the factors a , , bn), then we obtain for ClL a polynomial of degree (n - k ) + (n - (1 + L + 1 - k)) = 2n - 1 - ( 1 + L) in the a i and flj. The expansion of d e t ( C ( P , Q)) has the form: (*)
det(C(P, a))=
E
sg(°)Co, o(o)Cl,oo)"'" c,-1,o(,-1),
where S,, is the group of all permutations of the integers 0 , . . . , n - 1. N o w choose an arbitrary summand of (*), i.e. fix o ~ S,. Furthermore, we drop sg(o), since it is always equal to + 1 or 1. To complete the proof we need only show that -
deg(c0,,,(0)cl,o(1) ... c,,_ 1 , o ( n - 1)) = n2Since the degree of cz~ is 2n - 1 - ( I + L ) , it follows that deg(C,,o(,)) = 2 n - 1 - 1 - o ( I ) , and hence n--1
deg(c0,o(0)---Cn_l.o(n_l))
=
~ (2n-- 1--1--o(1)) I=0
*
n--1
=n(2n-1)-
n--1
E I-
Z o(1).
I=0
I=0
But Y'.z=oO(1) n-1 n--1 does not depend on o and is equal to Ei=0I. Hence, n--1
deg(Co,o
c._l,o(._l))=n(2n
-
1)-2
I I=0
----2n2-- n-- 2( n ( n -=- n1 2) ")2
[]
Lemma 6. C(P, Q) is singular if and only if P and Q have a common root. Proof. Assume x 0 is a common root. Then from the derivation of C(P, Q), we have the n × n system of homogeneous equations: n--1
~., CILX~=O,
O<.I<~n--1.
L=O
Since x ° = 1, we have a nontrivial solution which implies C(P, Q) is singular. Now assume
319
Y. De Montaudouirg tE. Tiller / The Cayley method
C(P, Q) is singular, i.e. det(C(P, Q ) ) = 0. Since a~ = flj implies det(C(P, Q ) ) = 0, it follows that ai - flj divides det(C(P, Q)). Furthermore, the ( a i - flj), 1 <~i, j <~n, are relatively prime in R [ a l , . . . , an, ill,-.., fin], and hence l-li~j=l(ai- flj) divides det(C(P, Q)). Since deg(I-li,j(ai flj)) = deg(det(c(P, Q)) = n 2, it follows from Lemma 5 that -
0 = det(C(P, Q))=pa~b~X--[(ai-flj), i,j where p is a nonzero integer constant. Therefore, a i = flj for some i, j.
[]
Proposition 7. I f P1,..., Pn are n polynomials which form a system of rank r in the basis 1, x .... , x ,, and if D is any nonzero polynomial, then DP1,.. ., DP n also forms a system of rank r. We now prove the fundamental theorem, Theorem 2: the rank loss of the Cayley matrix is equal to the degree of the gcd. Proof. Let P, Q ~ R[x], d e g ( P ) = d e g ( Q ) = n. Let D ~ R[x] be the gcd of P and Q, and deg(D) = r >I 0. Then P = DP 1, Q = DQ1, and deg(P 0 = deg(Q 0 = n - r. Using the above notation (derivation of the Cayley matrix), we see that D(xo)Pl(Xo)D(x)al(x ) - D(x)Pl(x)D(xo)al(xo) gp,e(x,
Xo) =
x - ~o
= D(x)D(xo)gP1,QI(X , Xo) n--r--1 n--r--1
=D(x)D(xo)
E
E
I=0 L=O
CILX~XI'
where [cIL ] is the Cayley matrix for P1, Q1- Lemma 6 implies that rank([CIL])=n--r. Multiplying through by D ( x ) and D(xo), we obtain
r(rx g,,,Q(~, Xo)-- ,~oE LE0=
)
ClLD(xo)~k D(x)~'.
So [C•L ] is also the Cayley matrix for P and Q in the basis D ( x ) , D ( x ) x .... , D ( x ) x n - ' - l . Rewrite gp,Q(X, XO) as gp,Q(X, XO)= [ D ( x ) , ' " , D ( X ) x n - r - 1 ] [ C I L ]
FID(x0). I
L D(xo)x~) - ' - 1 From Proposition 7 it follows that [ D ( x ) , . . . , D ( x ) x "-'-1] is independent and can be extended to a basis of the polynomials of degree n - 1:
~,Q(x, Xo)--1
fx(x0) = [1,/l(x),...,/,_l(x), D(x), D(x)x,...,D(x)x'-r-q[ei,]
fr-l(XO) D(Xo)
D(xo)xU "-1
Y. De Montaudouin,W. Tiller/ The Cayleymethod
320
where [elL] is the extension of [cIL] obtained by adding r null rows and columns. Clearly, rank([Si~ ]) = rank([cIL]) = n -- r. By a change of basis, we can find an n × n lower triangular invertible matrix A such that 1
D(x)
=
n--1
D(x)xn-,-] and it follows that
gp,Q(X, Xo)=[1, x, .... Xn-1]AT[cqL]AIlXo[ .
Hence, AT[a~L]A is the Cayley matrix for P and Q, and the observation that rank(AT[a~L]A) = rank([a~,]) = n - r completes the proof. []
The following corollary will be of fundamental importance in extracting the birational correspondence from the Cayley matrix (next section). We do not know if this corollary exists in the literature. Coronary 8. I f rank(C(P, Q ) ) = n - r , corner of C ( P , Q) is nonsingular.
then the ( n - r ) - t h
order submatrix in the lower right
Proof. I n t h e proof of Theorem 2 above we showed that C ( P , Q ) = A T C ( P 1, Q1)A, where A is lower triangular and invertible, and C(P1, Q1) is the extension of C(P1, Q1) which contains C(P1, Q1) in its lower right corner and zeros otherwise. Since P1 and Qt are relatively prime, it follows that C ( P 1, Qx) is nonsingular, and hence r a n k ( C ( P 1, Q 1 ) ) = rank(C(P1, Q1))= n - r. The following observation from linear algebra will complete the proof: let B be an arbitrary n × n matrix, and denote by Bp the (n - p ) - t h order submatrix in the lower right corner of B. Let C = ATBA, where A is lower triangular and invertible. Then for every p, 1 ~ p ~< n - 1, rank(Cp) = rank(Bp). But this can be seen by noting that, for i , j > p ,
Cij =
i.e. Cp =
E akibklalj = E akibklalj; l <~k,l<.n p+ l <~k,l<~n (aT)pBpap. []
Throughout this section we have assumed d e g ( P ) = deg(Q). We now show how to remove this restriction. Let P, Q ~ R[x], P ( x ) = •x"m " J, with m > n. Let a l , . . . , a m -,i=o a ix i, Q ( x ) = ~'~j=objx be the roots of P, and ill,---, fin the roots of Q. There are two ways to increase deg(Q) from n to m: (a) Regard Q as a polynomial of degree m, with the m - n leading coefficients equal to zero. (b) Multiply a by x " - n . Definition. The resultant of P and Q, denoted by Res(P, Q), is given by • ,.t \
re+n--1
rct n
m
re
amb il=il j----1 H
Y. De Montaudouin. W. Tiller / The Cayley method
321
Note that this definition does not require deg(P)= deg(Q). We will drop ( - 1 ) m+'-l, since it is not relevant to the discussion. It is easy to verify that
n (**)
ra
Res(P, Q ) = b ~ I - [ P ( f l / ) = a ~ , I " I a ( a i ) . j=l
i=l
When eliminating x, Res(P, Q) is actually what we are looking for when we set up the Cayley matrix; and if deg(P)= deg(Q), it is what we get, i.e. det(C(P, Q))= Res(P, Q) (see the proof of Lemma 6). If the degrees are not equal, we use either method (a) or (b) above. These methods yield a matrix C(P, Q) such that det(C(P, Q ) = F . Res(P, Q), F ~ R; that is, we get an extraneous factor. We now derive this factor. Definition. The reciprocal of a polynomial P, denoted by Rec(P), is the polynomial resulting
from reversing the order of the coefficients. If P ( x ) = a o + a m+
...
. . . +am xm, then R e c ( P ) =
+ a o Xm.
Proposition 9. I f ai are the roots of P, then 1 / a i are the roots of Rec(P). Lemma 10. Res(P, Q ) = Res(Rec(P), Rec(Q)). Proof. m
n
Res(Rec( P ), Rec(Q)) = agb~' Y1 I--I ( 1 / a i - 2/flj )
i=1 j=l
= agbg I-m[ f i (
i=1 j~l
flJ - ai ) Oli~j m
n
= ( a ~ b r ) ( a ~ n b r / a ; b r ) I-I I-I ( ¢ j - - a,) i=l
= Res(P, Q).
j=l
[]
Lennna 11. Let P, Q, T ~ R[x]. Then Res(P, QT) = Res(P, Q) Res(P, T). Proof. Follows from (* *).
[]
Now suppose we use method (a). Let Q denote the polynomial obtained by considering Q as a polynomial of degree m. Theorem 12. det(C(P, Q))= a,~-" Res(P, Q).
Proof. Clearly, R e c ( Q ) = x ~ - " Rec(Q). Then det(C(P, Q)) = Res(P, Q) = Res(Rec(P), Rec(Q)) = Res(Rec( P ), x m-" Rec( a )) = Res(Rec( P ), x ~ - " ) Res(Rec( P ), Rec( Q )) (from Lemma 11). From (* *) it follows that Res(Rec(P), x m - ' ) = (Rec(P(O))) m-" = a m-", and from Lemma 10, Res(Rec(P), R e c ( Q ) ) - R e s ( P , Q). Therefore, det(C(P, Q ) ) = amm - - n Res(P, Q). [] Now method (b). Let Q -- x m-'Q.
Y. De Montaudouin, W. Tiller / The Cayley method
322
T h e o r e m 13. d e t ( C ( P , Q ) ) = a on l - - r l Res(P, Q). Proof.
d e t ( C ( P , Q ) ) = R e s ( P , Q) = R e s ( P ,
xm-nQ)
= Res(P, x m-") Res(P, Q)= = a~'-" R e s ( P , Q).
(P(O))m-"Res(P, Q)
[]
So the extraneous factor, F, is just a~ ' - n if we use m e t h o d (b), a n d a,~ - " if we use m e t h o d (a). In a computer implementation, it should be checked which of the ai, i = 0, m, is simpler. T h e n d e t ( C ( P , Q)) must be divided by a 2 - " Again, we turn to Example 2 (see above). Consider P and Q as polynomials in x; i.e.
P, Q ~ C[yl[x]: 2
P(y, x ) = ( 7 y - 6 y 2 - y 3 ) + ( - 1
- 5y+ 3y2)x + x 2 =
E
ak x~,
k=0 1
a(y,
x)=(2-y+
10y2+y3)+(-2+
3 y - 7 y 2 ) x = ~_,btx'. l=0
m = d e g ( P ) = 2, n = deg(Q) = 1. We eliminate x. T h e 2 × 2 Cayley matrix has the form:
aobl - albo aob2- a2bo] aob2 - a2bo alb2- a2b l j" M e t h o d (a):
b 0 = 2 - y + 10y 2 + y 3 , bl = - 2 + 3y - 7y 2, b2 = 0.
D e n o t i n g the Cayley matrix by M ( y ) , we get:
M(Y)= [2_-25yy+_ 32y21 0 y 2 _ yX5y3 3 +
-2 +y - 10y 2 _y3]
2- 3y+7y 2
]
and d e t ( M ( y ) ) = - 1 2 y + 52y 2 - I45y 3 + 211y 4 - 205y 5 + 169y 6 - 70y 7. m--n = a~ = 1, it follows that d e t ( M ( y ) ) - R e s ( P , Q). Since a m M e t h o d (b):
bo = 0 b 1-- 2 - y + b2= -2+
10y 2 + y 3 , 3 y - 7 y 2.
D e n o t i n g the Cayley matrix by M ( y ) , we get: _~t(y) ---- [
14y -- 19y 2 + 78y 3 -- 54y 4 + 4y 5 -I-y 6
[ -14y
+ 33y 2 - 69y 3 + 45y 4 - 7y 5
-- 14y + 33y 2 -- 69y ~ + 45y 4 - 7y 5 ] 8y - 24y 2 + 43y 3 -- 21y 4
J
Y. De Montaudouin, W. Tiller / The Cayley method
323
d e t ( M ( y ) ) = _ 84y 2 + 436y 3 - 1339y 4 + 2399y 5 - 2846y 6 + 2624y 7 - 1709y 8 + 589y 9 - 70y 10 = (7y - 6y 2 + y 3 ) - d e t ( M ( y ) ) = a 0 • Res(P, Q).
3. T h e birational c o r r e s p o n d e n c e The birational correspondence is often stated a n d / o r motivated using rational curves in the xy-plane. In this context it can be stated as follows: let C(t)= (x(t), y(t)) be a parametric curve, where x(t), y ( t ) ~ K ( t ) , K a field. Assume that for each point on the curve there corresponds only o n e value of t. T h e n there exists f, g ~ K [x, y] such that t f(x, y)/g(x, y). In particular, if K is the field of rational numbers, then the coordinates (X(to), y(to)) of a point on the curve are-rational if and only if to is rational. In the next section we will show that this is a simple consequence of the Euclidean Algorithm; in this section we show that it is a result of eliminating t in the implicitization process (see implicitization example in the Introduction), a n d in particular, we derive from the Cayley matfix a formula which gives t as a function of x a n d y. The unit circle centered at the origin is the standard example for the birational correspondence. It is given by: (1)
x(t)=(1-t2)/(1
+ t 2)
y(t)=2t/(1 +t2).
and
Using substitution it is easy to verify that t = y / ( 1 + x). Let us derive this equation from the Cayley matrix. First, we rewrite (1) as:
F(x, y, t ) = ( x -
1)+(x+
1)t2=0,
G(x, y, t) = y - 2t + yt 2 = O. T h e n F, G E K[x, y][t]. Eliminating t, the Cayley matrix is: M(x, y ) = [ ( x - 1 ) ( - 2 ) - ( O ) ( y )
(x-1)(y)-(x+l)(y)]
(x-
1)(y) - (x + 1)(y)
-2y
2x + 2
__
]
( 0 ) ( y ) - (x + 1 ) ( - 2 )
E x p a n d i n g d e t ( M ( x , y)) and setting it equal to zero yields the circle. A n d f r o m -2y
2x + 2
x 2 + y2"= 1, the implicit-equation of
= 0,
it follows that ( - 2 y ) + (2x + 2)t --- 0, and hence, t = 2y/(2x + 2) = y / ( x + 1). N o t e that we could also have derived t = (1 - x ) / y (from the first line of M(x, y)). But at nonsingularities, these are equivalent since
y/(1 + x ) = y 2 / ( y + xy)= ( 1 - x 2 ) / ( y ( 1 = ((1 - x ) ( 1 +
x))/(y(1
+ x))=
+ x)) (1 -
x)/y.
T h e significance of the birational correspondence for the implicitization process is obvious. Suppose we have a parametric curve C(t) and a set of points {(x i, yi)}, 1 ~< i ~ n, and assume we want to k n o w if they lie on C(t), and if so, what are their corresponding t-values. This is a
Y. De Montaudouin, W. Tiller / The Cayley method
324
time consuming process when done directly from the parametric representation. Using Cayley we get the implicit form P(x, y ) = 0 of C(t) and the corresponding t = f ( x , y)/g(x, y). We substitute the (x;, y~) into P(x, y) to check if they lie on C(t), and if so, t = f / g gives their corresponding parameter values. N o w we look at the more general case, where the Cayley matrix may be larger than 2 x 2. For simplicity, we choose P, Q E K[x][y], K the field of rationals. The generalization to K[x a. . . . , x,][y] is obvious. Suppose d e g y ( P ) = degy(Q) = m > 0. Then
Lcm_l,0(x)
---
C..--1,m--,(X)
..--a
N o w suppose x0 is a value such that det[c~j(Xo)]= O. Then by Theorem 1, P(xo, y) and Q(x o, y), viewed as polynomials in K[y], have a gcd of degree = r > 0, where r is the rank loss of [cij(xo)]. Suppose r = 1. By Corollary 8, we let M 1 denote the ( m - 1)-th order nonsingular submatrix in the lower right corner of [cij(Xo)]. Then
Fy0] 1 c ,o(xo)
(2)
~
"
o
~
~
Lyg,-
.
LCm-l,o(Xo)
and we can solve for Y0 using Cramer's rule: det(A) Yo = d e t ( M 1 ) ' where A is the matrix formed by replacing the first column of M 1 with the right hand side of (2). Clearly, Y0 is rational if x0 is rational. Obviously, we cannot use the above method to solve for Yo if r > 1. In this case there is no birational correspondence. For curve implicitization r = deg(gcd) > 1 implies a singular point on the curve; for curve/curve intersection it implies either that the curves touch tangent to one another or that there are several y-solutions corresponding to one x-solution. We now return to Example 2 of the previous section. W h e n we substitute x o = 1 into the Cayley matrix we obtain: M(1) =
[i 0 0] 12 0
0 • -6
Then r a n k ( M ( 1 ) ) = 2. The last two columns are independent, and if we discard the first row ( k = 1), we get: M I = [12 0
O] -6 "
[1~
Y0 =
Hence, 0
6][Yogi [:] _
and Y°=-
I°0
-6
0
° l=0
-6
Y. De Montaudouin, W. Tiller / The Cayley method
325
This shows that (1, 0) is a solution. For x 0 = 2, rank(M(2))= 1, which implies there is no birational correspondence in this case. This is clear since the gcd of P(2, y) and Q(2, y) is quadratic, i.e. it could have complex roots.
4. The connection between the Cayley method and the Euclidean Algorithm
Let t'1, P2 ~ C[x, y]. We may consider/'1 and P2 as belonging to C(x}[y]; i.e. they are polynomials in the variable y with coefficients which are rational functions in x. N o w C( x } is a commutative field, which implies that C(x }[y] is a Euclidean domain. It follows that we have a Euclidean algorithm in C( x }[y]. Furthermore, C( x }[y] a Euclidean domain implies that it is a principal domain, and hence, we have the following result: if P1 and P2 have no common components (factors) in C[x, y], then they are relatively prime in C(x }[y] and there exist U'1, U~ ~ C{x}[y] such that: U~P1+ U~P2= 1. Now if R is a common multiple of the denominators of the coefficients of U{ and U~, then we can multiply through by R and obtain: U1P1 + U2P2= R, where R ~ C[x] and U1, U2, P1, P2 ~ C[x, y]. Now to find R we can do a Euclidean algorithm in C(x}[y]. Assume degy(P2)>/degy(P1) > 0. Dividing P2 by P1 we get: P2 = P1Q1 + R1, Q1, R1 ~ C(x}[y] and degy(R1) < degy(P1). If we repeat this enough times, we eventually get a remainder equal to zero. For illustration, suppose this requires four steps. Then: P2 = P I Q 1 + R1,
P1 = RaQ2 + R2, R 1 = R 2 Q 3 + R 3, R 2 = R 3 Q 4 + O.
Now R 3 is a constant in C{x)[y], that is R 3 ~ C ( x ) , and we have: R 3 = R 1 - R2Q 3
=Ra-Q3(P1-R1Q2 ) = R1 -
Q3PI+
Q302(P2 - P1Qa)
= R 1 -- Q 3 P 1 + Q 3 Q 2 P 2 - P I Q I Q 2 Q 3 = P2 - P I Q 1 - Q 3 P 1 + Q 3 Q 2 P 2 - P 1 Q I Q 2 Q 3
= P I ( - - Q 1 - Q3 - Q 1 Q 2 Q 3 ) + / ' 2 ( 1 + a 2 0 3 ) ,
where Qa, Q2, Q3 ~ c{x}[y] and R 3 ~ C ( x ). Now multiplying through by a common multiple of the denominators of R3, 01, Q2 and Q3 yields: R=P1U1 +P2U2, where U1,
U2 ~ C[x, y] and R ~ C[x]. When we eliminate y from Pa(x, y) and P2(x, y), by the Cayley or any other method, it is R that we find. But R2 is also of interest, as it gives us an expression of y as a function of x. Since d e g y ( R 3 ) = 0, it follows that dege(Rz)>~ 1. Suppose d e g y ( R 2 ) = 1. Then R2(x, y)= Ao(x)+yAa(x), where Ao, A I ~ C ( x } . N o w if Xo is a root of R, and hence, of R3, then R2(xo, y) is the gcd of Pl(xo, y) and P2(x o, y). Hence, if Yo is a root of R2(xo, Y), then (Xo, Yo) is a solution to the system Pl(X, y)=P2(x, y ) = 0 , and Rz(xo, yo)=Ao(xo)+ yAa(xo)--- 0, which implies that y = --Ao(xo)/Aa(xo). This expression is the birational correspondence, and it is equivalent to the expression of y as a ratio of two sub-determinants derived above as a by-product of the Cayley matrix, in the sense that they yield the same value of y for those values of x which are x-solutions to the system Pl(x, y)= P2(x, y ) = 0. The equation y = - A J A x corresponds to a rank loss of one in the Cayley matrix. If there are several
326
Y. De Montaudouin, W. Tiller / The Cayley method
y-solutions for some x-solution, then the gcd R 2 ( x , y) can be of degree k > 1 in y: R z ( x , y ) = A o ( x ) + .. • + y k A k ( x ). This corresponds to a rank loss of k in the Cayley matrix.
Acknowledgement We thank Ron Goldman for reading the original manuscript and making several valuable suggestions, one of which led to the current, stronger version of Corollary 8.
References Griffiths, H.B. (1981), Cayley's version of the resultant of two polynomials, Amer. Math. Monthly 88 (5), 328-338. Jacobson, N. (1974), Basic Algebra 1, Freeman, San Francisco. Kajiya, J.T. (1982), Ray tracing parametric patches, Computer Graphics 16 (3), 245-253. Lelong-Ferrand, J. and Arnaudies, J.M. (1978), Cours de Mathematiques, Tome 1 Algbbre, Dunod, Paris. Sederberg, T., Anderson, D. and Goldman, R. (1984), Implicit representation of parametric curves and surfaces, Computer Vision, Graphics and Image Processing 28, 72-84. Sederberg, T., Anderson, D. and Goldman, R. (1985), Implicitization, inversion, and intersection of rational cubic curves, to appear. Uspensky, J.V. (1948), Theory of Equations, McGraw-Hill, New York.