This section describes an algorithm, which uses generalized eigenvalues and eigenvectors, for solving a system of two bivariate polynomials , with . It treats multiple roots, which are usually considered as obstacles for numerical linear algebra techniques. The numerical difficulties are handled with the help of singular value decomposition (SVD) of the matrix of eigenvectors associated to each eigenvalue. The extension to equations could be treated similarly.
A vector is called a {generalized eigenvector} associated to the eigenvalue if . The matrices and have generalized eigenvalues if and only if . If , then can be finite, empty, or infinite. Note that if then . Moreover, if is invertible then , which is the ordinary spectrum of .\
Recall that an matrix with polynomial coefficients can be equivalently written as a polynomial with matrix coefficients. If , we obtain , where are matrices.
Definition: The companion matrices of are the two matrices and given by:

And we have the following interesting property :
Proposition: With the above notation, the following equivalence holds:

Hereafter we will consider systems having only a {finite} number of roots. This condition is not quite restrictive since it is sufficient (and necessary) to require that polynomials and are coprime in (otherwise one can divide them by their ).
In the case where one of the curve is a line, i.e. one of the polynomial has degree 1, the problem is reduced to solving a univariate polynomial. Indeed, one may assume e.g. that and thus we are looking for the roots of . Let us denote them by . Then we know that
where the 's are assumed to be distinct and is a nonzero constant in . The integer , for all , is called the {multiplicity} of the root , and it turns out that if do not vanish at infinity in the direction of the axis (in other words, if the homogeneous part of of highest degree do not vanish when ). This later condition can be avoid in the projective setting: let be the homogeneous polynomial obtained by homogenizing with the new variable , then we have
where is an integer corresponding to the multiplicity of the root at infinity, and . Moreover, it turns out that the roots and their corresponding multiplicities can be computed by eigenvalues and eigenvectors computation.
In the sequel we will generalize this approach to the case where and are bivariate polynomials of arbitrary degree. For this we first need to recall the notion of multiplicity in this context. Then we will show how to recover the roots from multiplication maps.
A modern definition of the intersection multiplicity of and at a point is (see [Ful84])
where the point is here abusively (but usually) identified with its corresponding prime ideal in , denoting the ordinary localization of the ring by this prime ideal. As a result, the finite algebra (which is actually finite if and only if and are coprime in ) can be decomposed as the direct sum
and consequently
The intersection multiplicities can be computed using a resultant. The main idea is to project ``algebraically'' the intersection points on the axis. To do this let us see both polynomials and in where , that is to say as univariate polynomials in with coefficients in the ring which is a domain; we can rewrite
Their resultant (with respect to ) is an element in which is nonzero, since and are assumed to be coprime in , and which can be factorized, assuming that and do not vanish simultaneously, as
where the 's are distinct elements in and (as sets). For instance, if all the are distinct then we have and for all . Moreover, we have (see e.g. [Ful84] (S)):
Proposition: For all , the integer equals the sum of all the intersection multiplicity of and at the points such that . 
As a corollary, if all the 's are distinct (this can be easily obtained by a linear change of coordinates ) then is nothing but the valuation of at . Another corollary of this result is the wellknow Bezout theorem for algebraic plane curves, which is better stated in the projective context. Let us denote by and the homogeneous polynomials in obtained from and by homogenization with the new variable .
Proposition: [Bezout theorem] If and are coprime then the algebraic projective plane curves associated to and intersect in points in , counted with multiplicities. 
{proof} It follows straightforwardly from the previous proposition and the wellknown fact that is a homogeneous polynomial in of degree (see e.g. [Lang02], [VWe48II]). {proof}
Suppose given such a point and assume that and are not both zero. If is of dimension one, then it is easy to check that there is only one point such that and moreover any element in is a multiple of the vector . Therefore to compute , we only need to compute a basis of and compute the quotient of its second coordinate by its first one. In case of multiple roots, this construction can be generalized as follows:
In this theorem, observe that the construction of and is always possible because we assumed that both polynomials and depend on the variable , which implies that the number of rows in the matrix is always strictly greater than the .
The previous theorem shows how to recover the coordinates of the roots of the system from the Sylvester matrix. However, instead of the Sylvester matrix we can use the Bezout matrix to compute (see section [*] ). This matrix have all the required properties to replace the Sylvester matrix in all the previous results, except in proposition [*] . Indeed, the Bezout matrix being smaller than the Sylvester matrix, the condition given in remark [*] is not always fulfilled; using Bezout matrix in proposition [*] requires that for all roots of ,
This condition may failed: take for instance in the system
Note that the use of the Bezout matrix gives, in practice, a faster method because it is a smaller matrix than the Sylvester matrix, even if its computation takes more time (the save in the eigencomputations is greater).
For difficulties ``c'', it is necessary to use the singular value decomposition: If a generalized eigenvalue of and , and the matrix of the associated eigenvectors, then the singular value decomposition of is
and is the rank of . Let
The matrix is also a matrix of eigenvectors, because the changes are done on the columns, and V is invertible. Hence and have the same rank , and is a matrix of linearly independent eigenvectors. Consequently can be used to build and .
This leads to the following complete algorithm:
{Algorithm { {Numerical approximation of the root of a square bivariate system} {{0.3cm}}
{Input:} Two bivariate polynomials
template<class C> Seq<VectDse<C> > solve(p,q,EigenRes<C>());
synaps/resultant/bivariate_res.h
.#include <synaps/msolve/projection/EigenRes.h> #ifdef SYNAPS_WITH_LAPACK using std::cout; using std::cin; using std::endl; typedef MPol<double> mpol_t; int main(int argc, char **argv) { mpol_t p1("16*x0^4160*x0^2*x1^2+400*x1^432*x0^3+160*x0*x1^2+6*x0^250*x1^2+10*x0+1.5625"), p2("x0^2x0*x1+x1^21.2*x00.0625"); cout<<"p1="<<p1<<endl; cout<<"p2="<<p2<<endl; Seq<VectDse<double> > res; Seq<int> mux,muy; res= solve(p1,p2, EigenRes<double>(1e3)); res= solve(p1,p2,mux, EigenRes<double>(1e6)); res= solve(p1,p2,mux,muy, EigenRes<double>(1e7)); cout<<"Intersection points:"<<endl; for(unsigned i=0;i<res.size();++i) std::cout<<"("<<res[i][0]<<","<<res[i][1]<<") of multiplicity "<<mux[i]<<" "<<muy[i]<<std::endl; } #else #include <iostream> int main() { std::cout <<" LAPACK is not connected"<<std::endl; return(77); } #endif //SYNAPS_WITH_LAPACK