\documentstyle[11pt]{article}\topmargin-.25in\oddsidemargin.05in\evensidemargin.05in\textheight8.55in\textwidth6in\newtheorem{theorem}{Theorem}[section]\newtheorem{proposition}{Proposition}[section]\newtheorem{lemma}{Lemma}[section]\newtheorem{corollary}{Corollary}[section]\newtheorem{definition}{Definition}[section]\newtheorem{algorithm}{Algorithm}[section]\newtheorem{remark}{Remark}[section]\renewcommand{\baselinestretch}{1.5}\renewcommand{\theequation}{\thesection.\arabic{equation}}\begin{document}\title{Computation of Approximate Polynomial Gcds and an Extension.}\author{   Victor~Y.~Pan \\Mathematics and Computer Science Department \\ Lehman College, City University of New York\\ Bronx, NY 10468\\Internet:  VPAN@LCVAX.LEHMAN.CUNY.EDU \\(Supported by NSF Grant CCR 9625344 \\and PSC CUNY Award No. 667340) } \maketitle\begin{abstract} Computation of approximate polynomial gcds  is important both theoretically and due to applications, inparticular to linear control systems. We study two approaches to the solution so far omitted by theresearchers, in spite of intensive recent work in this area. Correlation to numerical Pad\'{e} approximation enabled us toimprove computations for both problems (gcds and Pad\'{e}). Reduction to approximating polynomial zeros enabled us to obtaina new insight into the gcd problem and to devise effectivesolution algorithms. In particular, unlike the known algorithms,we compute an upper bound on the degree of approximate gcds at a lowcomputational cost, and this enables us to certify the correctness ofthe solution. We also argue in favor of restating the problem interms of the perturbation of the zeros (rather than coefficients) of theinput polynomials,  which leads us to the solution via the computation ofmaximum matchings or connected components in the associated bipartitegraphs. \vspace{.15in}\noindent{\bf Key words:} polynomial gcds,approximate gcds, Pad\'{e} approximation, Hankelmatrices, polynomial zeros, root neighborhoods, bipartite graphs. \vspace{.15in}\noindent{\bf 1991 Mathematics Subject Classification:} 68Q40, 65D99, 65Y20.\end{abstract} \section{Introduction.}\label{sec1} \setcounter{equation}{0}  Computation of polynomial {\it greatest common divisors (gcds)} is afundamental problem of algebraic computing and has important widespread applications in network theory and control linear systems (see \cite{b81}, \cite{k80}, \cite{k87}, and [\cite{pb73}), which require numerical solution of the problem, where the input is givenapproximately, within some fixed error bounds.  On the other hand,computation of polynomial gcds is  an excellent example of numericallyill-posed problems.  For instance, let $v(x)$ be a non-constant divisor of a polynomial $u(x)$.   Then gcd $(u(x),v(x))=v(x)$, but gcd $(u(x)+\delta ,v(x))=1$ for any constant $\delta\neq 0$. Thus, a small perturbation of$u(x)$ may cause a {\it dramatic decrease} of the degree of the gcd.\cite{sc85}, \cite{ns91}, \cite{km94},\cite{hs95}, \cite{cgtw95}, \cite{kl96}, and \cite{egl97} define  approximate gcds soas to avoid the latter deficiency. Namely, for two polynomials$u(x)$ and $v(x)$, a fixed polynomial norm, and a positive $\epsilon$,  one may non-uniquely define an approximate gcd or, we say {\it $\epsilon$-gcd}  $d^* (x)$,  as gcd $(u^* (x),v^* (x))$ whose degree $d^*$ is maximized over allpolynomials$u^*(x)$ and $v^*(x)$ in an $\epsilon$-neighborhood of $u(x)$ and $v(x)$, satisfying \begin{equation}\label{equ1.1} \ {\rm deg} \  u^* (x)\leq \ m = \ {\rm deg} \   u(x) , \ \ {\rm deg} \   v^* (x)\leq \ n =  \ {\rm deg} \   v(x) , \end{equation}\begin{equation}\label{equ1.2} \parallel u^* (x)-u(x)\parallel \leq  \epsilon   \parallel u(x) \parallel , \ \ \\parallel v^* (x)-v(x)\parallel \leq  \epsilon  \parallelv(x)\parallel .\end{equation} In spite of extensive work and substantial progress, there remain several open problems with the computation of $\epsilon$-gcds. The pioneering paper \cite{sc85} studies only the asymptoticcomplexity of computing the $\epsilon$-gcds where $\epsilon\rightarrow 0$. The results of \cite{sc85} rely on extremely tedious analysis but donot apply to the most realistic and practically important case, where theinput errors can be relatively large,  and the precision of computing issufficiently high to ignore rounding errors of the computation. The latterassumptions  imply the model of study with inexact input and infinite precision computations. In our paper, we will assume thismodel, following all  papers cited above except for \cite{sc85}. On the other hand, all these papers, except for \cite{kl96}, presentheuristic solution algorithms, whereas the arithmetic complexityestimates for the solution presented in \cite{kl96} are quite large(that is, unspecified polynomial in $m+n$ and exponential in $d^*$), thus implying the need for further study. Furthermore, neither of theso far presented algorithms allows its effective parallelization. Technically, there is another major omission. Several approaches havebeen studied or at least listed so far. They rely on Euclideanalgorithm  \cite{sc85}, \cite{hs95}, \cite{egl97}, computations with subresultant matrices \cite{cgtw95}, \cite{egl97}, various techniques of least-squares computations,optimization, and quadratic programming \cite{cgtw95}, \cite{kl96}, and a version of Lazard's algorithm [L81], which is equivalent to the matrix pencil algorithm of [KM94] (described in [KM94] by using the terminology of automatic control). At least two important approaches are missing, however, from all these papers, that is, ones based on {\it Pad\'{e} approximation} and {\it approximating polynomial zeros}. In our present study, wedemonstrate several methodological and computational advantages of thesetwo approaches over the cited ones.  Pad\'{e} approximation is an important and well developed subjectwith applications to algebraic computing, signal processing, and thestudy of analytic functions \cite{gr72}, \cite{bgy80}, \cite{bp94}. We show some computational benefits of reducing theapproximate gcd problem to the computation of Pad\'{e} approximations for an inexact input, versus the subresultant approach, intensivelystudied in  \cite{cgtw95}, \cite{kl96}. (In particular, the formerapproach involves better structured matrices of smaller size.) On the other hand, the reduction into the opposit direction, from Pad\'{e} computations to the gcd computations, enables us tosolve both problems by our second approach, which, as weshow, is quite effective. Namely, we reduce the gcd problem to approximating the zerosof the input polynomials, where highly effective algorithmsare available \cite{p95}, \cite{p96}, \cite{p97}. With such a reduction, we explain the numerically unstablebehavior of approximate gcds and restate the problem in terms of theconcept of the {\it polynomial root neighborhoods} of  \cite{m86}. This enables us to choose a proper class of theinput perturbations for the gcds and also to associate {\itappropriate bipartite graphs} with the gcds. Then the solution of theapproximate gcd problem is reduced to computing {\it maximummatchings} in such graphs. In particular this enables us to computeeasily a non-trivial upper bound on the degree of the approximate gcdsand to certify correctness of the solution, which was the bottleneck ofthe known approaches. {\it We made this progress assuming the customary formulation of thegcd problem}, based on the {\it perturbation of the inputcoefficients}. In sections \ref{sec9} and \ref{sec10} we study theproblem based on the {\it perturbation of the zeros} of the inputpolynomials and argue that this is a more appropriate basicassumption. We also show its computational advantages: it alwaysenables us to compute an approximate  gcd itself (rather than acandidate polynomial, which may have a lower or higher degree) andto simplify the stage of the computations in bipartite graphs, byreplacing the matching stage by the computation of the {\itconnected components}. The proposed algorithms have lower computational complexity (which rangesbetween linear and quadratic in terms of arithmetic operations andcomparisons involved, except for$O(n^{2.5})$ comparisons used for matching), and their parallelization enables NC (or RNC) and work efficient solution of theapproximate gcd problem (cf. e.g. \cite{j92} or \cite{bp94}, ch. 4, onthe definitions of NC, RNC and work efficiency). We organize our presentation as follows. After some preliminatires insection \ref{sec2}, we study the Pad\'{e} approximation approach insections\ref{sec3} and \ref{sec4}. In sections \ref{sec5}--\ref{sec7} westudy the $\epsilon$-gcds problem based on the concept ofa polynomial root neighborhood. In sections \ref{sec8} and\ref{sec9}, we study the$\delta$-gcd problem,  where the perturbation is applied to the zerosof the input polynomials. In section \ref{sec10}  we comment on theextension of our approach to Pad\'{e} approximation for an inexactinput. Section\ref{sec11} is left for a short discussion. In the appendix we brieflyrecall some major known methods for $\epsilon$-gcds and prove anauxiliary result for testing $\epsilon$-divisibility. \vspace{.52in}{\bf Acknowledgements.} I gratefully acknowledgereceiving reprints of [CGTW95]from Andr\'{e} Galligo and Erich Kaltofen, a preprint of [HS95] fromthe former, and a reprint of [HZ96] from Uri Zwick.My present work was substantially motivated by  Andr\'{e} Galligo's comments on [CGTW95], delivered in [G95], and partly by Erich Kaltofen's interest to thecomputaton of numerical rank of a Toeplitz matrix. This paper was written and submitted for publication in 1995, turned intoa research report\cite{p96a} in 1996, and substantially revised in 1997, based on thevaluable suggestions by the referees.\vspace{.2in}\section{Polynomial and vector norms,$\epsilon$-divisibility and $\epsilon$-gcds (some definitions and preliminaries).}\label{sec2}  \setcounter{equation}{0}  Hereafter, we will refer to arithmetic operations as to {\it ops}.A polynomial $p(x) = \sum_{i=0}^k p_i x^i$ can be identified with its coefficient vector $\vec{p} = [p_0, \ldots, p_k]^T$. (We write $w^T$ for the transpose of a vector or a matrix $w$.) The same norms will be used for both vectors and polynomials, in particular to measure the distances, dist$(s(x), t(x))= \| s(x)-t(x)\|= \| \vec{s}-\vec{t}\|$. We recall the customary  norms \begin{equation}\label{equ2.1}\| \vec{p} \|_h = \| \sum_i p_i x^i \|_h = (\sum_i|p_i|^h)^{1/h}, \end{equation} which for $h=\infty$ turn into the {\it maximum norm}, $$ \| \vec{p} \|_{\infty} = \| \sum_i p_i x^i \|_{\infty} =\max_i|p_i|. $$ In some cases, the weighted $h$-norms, $\| \vec{p} \|_{h, \vec{w}} = \| \sum_i p_i x^i \|_{h, \vec{w}} =(\sum_i|p_i|^h w_i)^{1/h}$ for a fixed vector $\vec{w}=(w_i)$, are also useful in the study of $\epsilon$-gcds \cite{egl97}.  When our study applies to any fixed norm, we will write$\| . \|$. The next two defintions re-introducethe $\epsilon$-gcds of (\ref{equ1.1}),(\ref{equ1.2}), based on the concept of an $\epsilon$-divisor. \begin{definition}\label{def2.1}  A polynomial$d(x)$ is an {\it$\epsilon$-divisor} of  a polynomial $p(x)$ (under afixed norm $\|.\|$) if there exists a perturbation of $p(x)$ by a polynomial $\Delta (x)$ such that $d(x)$ divides $p(x) +\Delta(x)$, $\deg \Delta(x) \leq N$, and $\| \Delta(x) \|\leq \epsilon \|p(x) \|$.\end{definition}\begin{definition}\label{def2.2} For two polynomials,\begin{equation}\label{equ2.2} u(x) = \sum_{i=0}^m u_i x^i, \  \   v(x) = \sum_{i=0}^n v_i x^i, \  \  u_m v_n \neq 0, \  \  m \leq n, \end{equation}every their monic $\epsilon$-divisor $g(x)$ that has the maximum degree, $d_{\epsilon}=d(u,v,\epsilon)$, is called their {\it $\epsilon$-gcd}. \end{definition}Using the {\it 2-norm}, $\| .\|_2$, often leads to some computational advantages. \begin{proposition}\label{pro2.1}\cite{cgtw95}.  Given two polynomials $p(x)$ and $d(x)$ of degrees $k$and $l$,  respectively, $k>l$, and  a positive $\epsilon$, it sufficesto use $O(k \log (k-l) + \min\{(k-l) \log^2(k-l), kl\})$ ops to decide  if $d(x)$ is an {\it $\epsilon$-divisor} of $p(x)$ (under the 2-norm). \end{proposition}The algorithm of \cite{cgtw95} supporting this proposition actuallycomputes two polynomials $q(x)$ and $\Delta (x)= p(x)-d(x) q(x)$ suchthat $\|\Delta (x)\|_2$ is minimum, and $\deg \Delta (x) \leq  \deg p(x)$. Having $\Delta(x)$ available, we may check immediately if  $\|\Delta(x)\|_h \leq \epsilon \| p(x) \|_h$ for$h=1$ or $h = \infty$; the values $\|\Delta(x)\|_h$ for $h=1$ arewithin factors $n$ or $\sqrt{n}$ from the minimum $\| p(x) -d(x) q(x)\|_h$,   due to the following useful relations: $$ \| \vec{v} \|_{\infty} \leq \| \vec{v} \|_2 \leq \| \vec{v} \|_1 \leq\sqrt{n} \| \vec{v} \|_2 \leq n \| \vec{v} \|_1, $$ which hold for any vector $\vec{v}$ of dimension $n$. Furthermore, we can prove the following extension ofproposition~\ref{pro2.1} (see appendix \ref{appb}). \begin{proposition}\label{pro2.2} The $\epsilon$-divisibility under aweighted  2-norm can be decided at the cost of performing $O(N^2 \log N)$ ops.\end{proposition} \begin{remark}\label{rem2.1} Our study can be easily extended to $(\epsilon_{\vec{u}}, \epsilon_{\vec{v}})$-gcds where $g(x)$ is an  $\epsilon_{\vec{u}}$-divisor for $u(x)$ and an$\epsilon_{\vec{v}}$-divisor for $v(x)$, for fixed $\epsilon_{\vec{u}}$and $\epsilon_{\vec{v}}$, or more generally, where a distinctupper bound may be imposed on the allowed perturbation of eachcoefficient of $u(x)$ and $v(x)$.\end{remark}   \section{Pad\'{e} approximation approach to computing a polynomial gcd or $\epsilon$-gcd.}\label{sec3}\setcounter{equation}{0}  A pair of polynomials $w(x)$ and $z(x)$ of the smallestdegree that satisfy the equation $w(x) v(x) = z(x) u(x) $ immediately defines the gcd, $g(x) = {\rm gcd} \ (v(x) , u(x)) $, as the polynomial $$ g(x) = \frac{u(x)}{w(x)} =  \frac{v(x)}{z(x)}. $$ The problem of the computation of $w(x)$ and $z(x)$ can   be appropriately re-phrazed by using the important classical concept of Pad\'{e} approximation (cf. \cite{gr72},  \cite{bgy80},\cite{bp94}). \begin{definition}\label{def3.1}  For any formal power series $a(x)=\sum_{i=0}^{\infty} a_i x^i$ and for two non-negative integers $k$ and $\ell$, a pair of polynomials $s(x)$ and $t(x)$ is a $(k,l)$-th {\it Pad\'{e} approximation} of $a(x)$ if $\deg s(x) \leq k$, $\deg t(x) \leq l$, and  $s(x) - a(x)t(x) = 0 \bmod x^{N+1}$, $N=k+l$. \end{definition}\begin{proposition}\label{pro3.1} \cite{gr72}. The pair of polynomials $q(x)$ and $t(x)$ of definition~\ref{def3.1} is defined uniquely, up to its scaling by common factors or common divisors.\end{proposition} Given the polynomials $u(x)$ and $v(x)$ of (\ref{equ2.2}), we mayassume that $v(0)\neq 0$ and define the formal power series \begin{equation}\label{equ3.1} h(x) =\sum_{i=0}^{\infty} h_i x^i = \frac{u(x)}{v(x)}.\end{equation} (The restriction $v(0)\neq 0$ can be removed by removing the maximum degree factors $x^i$ and $x^j$ from $u(x)$ and $v(x)$, respectively, andadding the factor $x^{\min\{i,j\}}$ to their gcd. Alternatively, one maywork with the reverse polynomials $U(x) = x^m u(1/x)$, $V(x) = x^n v(1/x)$, $W(x) = x^{m-d}w(1/x)$, $Z(x)= x^{n-d} z(1/x)$.) The computation of $h(x) \bmod x^{N+1}$ amounts to the computationmodulo $x^N$ of a  polynomial reciprocal and a polynomial product andcosts $O(n \log n)$ ops, \cite{bp94}, p. 22. The $(k, l)$-th Pad\'{e}approximation can be computed by means of the extended Euclideanalgorithm at the cost $O(N \log^2 N)$ ops for $N=k+l$, \cite{bgy80}, \cite{bp94}, pp. 38--39. These considerations give us a gcd algorithm, and we may extend it to computing an $\epsilon$-gcd as follows. \begin{description} \item[a)] For each $d$, $d=0, 1, \ldots, m$, successively orconcurrently compute $w_d(x)$, $z_d(x)$, the $(m-d, n-d)$-th  Pad\'{e}approximation to $h(x)$; then apply proposition \ref{pro2.1} to test if$w_d(x)$ is an $\epsilon$-divisor of $u(x)$. If "not", discard sucha value $d$. Otherwise compute a  polynomial $g_{d, \epsilon}(x)$ (an $\epsilon$-quotient) such that $\| w_d(x) g_{d,\epsilon}(x) -u(x)\| \leq \epsilon$.  Then applyproposition \ref{pro2.1} again, to test if $g_{d, \epsilon}(x)$  is an$\epsilon$-divisor of$v(x)$ too. If "not", discard this $d$.  Otherwise, store such a $d$ and $g_{d, \epsilon}(x)$ . \item[b)] If  the set of $d$'s is empty, output FAILURE; otherwiseselect the largest remaining $d=d^*$ and output $d^*$ and $g^*(x) =g_{d^*, \epsilon}(x)$ . \end{description} By construction, $g_{d^*, \epsilon}(x)$ is a common $\epsilon$-divisor of$u(x)$ and $v(x)$, though some $\epsilon$-perturbations of $u(x)$ and $v(x)$ may have a common $\epsilon$-divisor of a larger degree. The algorithm is performed at quite a low computational cost of $O(mn \log^2 n )$ ops. Furthermore, the stage of computing the $(m-d, n-d)$ Pad\'{e} approximations of $h(x)$ for all $d$ requires only $O(n \log^2 n)$ ops, except for some cases of degeneration\cite{bgy80}, sect.~5. On the other hand, this approach, like the Euclidean and subresultant approaches, gives us no effective general means for detecting degeneration or verifyingthat the output polynomial is an $\epsilon$-gcd. \section{The Hankel/Bezout techniques for Pad\'{e}approximation.}\label{sec4} \setcounter{equation}{0}  In the algorithms of the previous sections, one may compute  Pad\'{e}approximations by relying on computations with Hankel matrices, rather than on the extended Euclidean algorithm. This approach may require  a little more ops but  enables better numerical control of the computations and provides an additional insight into the behavior of $\epsilon$-gcds. The basic idea is to associate   the formal power series $h(x)$ of (\ref{equ3.1}) with the infinite Hankel matrix $H= H(u,v) =(h_{i,j})$, $ h_{i,j}= h_{i+j+m-n+1}$,  $i, j =0, 1, \ldots$. For readers' convenience, we recall that a matrix $(h_{i,j})$ is called a {\it Hankel matrix} if all its entries are invariant intheir shift into antidiagonal direction, that is, if  $ h_{i,j}= h_{i+1, j-1}$ for all pairs $i,j$ for which $ h_{i,j}$ and $ h_{i+1, j-1}$  are defined. A general $n \times n$ matrix may have $n^2$ distinct entries, but an $n \times n$ Hankel matrix is symmetric  ($ h_{i,j}= h_{j,i}$), has at most $2n -1$ distinct entries, and is completely defined by the pair of its first and last rows orcolumns. Computations with Hankel matrices are also dramatically simplifiedversus the case of general matrices. In particular we have the following results \cite{bp94}. \begin{proposition}\label{pro4.1} The multiplication of an $N \times N$ Hankel matrix by a vectorcan be reduced to multiplication of two polynomials of degree $O(N)$and performed in $O(N\log N)$ ops. \end{proposition} \begin{proposition}\label{pro4.2}$O(N\log^2 N)$ ops suffice to solve a non-singular linearsystem of $N$ equations with a Hankel coefficient matrix.  \end{proposition} Now let $H_k$ denote the $k \times k$ leading principal submatrix  of $H$, that is, one where $i$ and $j$ range from 0 to$k-1$.  (Clearly, $H_k$ is a Hankel matrix.) We have the following results (\cite{gr72}, \cite{bgy80},\cite{bp94}, pp. 140). \begin{proposition}\label{pro4.3} Let $r = {\rm rank} \  H_n$, $d= {\rm gcd} \  g(x)$. Then the matrix$H_r$ is non-singular, and $d+r =n$. \end{proposition} \begin{proposition}\label{pro4.4} The computation of the $(m,n)$ Pad\'{e} approximationof the formal power series  $h(x)$ of (\ref{equ3.1})can be reduced to the solution of a consistent linearsystem of equations with the coefficient matrix$H_n$, multiplication of an $m \times m$triangular Hankel matrix by a vector, and asubstraction of a pair of $m$-dimensional vectors from each other. \end{proposition} The harder part of the latter computation is the solution of the linear system. It can be performed in three steps: \begin{description} \item[a)] compute the rank $r$ of the matrix $H_n$, \item[b)] solve a non-singular linear system of equations with thecoefficient matrix $H_r$ (non-singularity is byproposition~\ref{pro4.3}, the solution cost is $O(r \log^2 r)$ byproposition~\ref{pro4.2}). \item[c)] recover the solution of the original consistent linear systemwith the matrix $H_n$ (in $O(n \log n)$ ops by proposition \ref{pro4.1}). \end{description} \begin{remark}\label{rem4.1} At stage b), one may substantially improve numerical stability of the computations by means of a simple transitionfrom the Hankel to Bezout linear systems (\cite{bp94}, p.162). \end{remark}  In the context of the computation of $\epsilon$-gcds, the matrix $H_n$should be allowed to vary with the input polynomials $u(x)$ and $v(x)$.The $\epsilon$-perturbations of $u(x)$ by $\delta_u(x)$ and  $v(x)$ by $\delta_v(x)$, $ \| \delta_u(x)\| \leq \epsilon$,$\|\delta_v(x)\leq \epsilon$, cause the perturbations of the polynomial $(u(x)/v(x)) \bmod x^{m+n+1}$  within $\delta_{u, \epsilon} + \delta_{v, \epsilon}$ where \begin{equation}\label{equ4.1} \delta_{u, \epsilon}= \epsilon \| (1/v(x)) \bmod x^{m+n+1}\|, \end{equation} \begin{equation}\label{equ4.2} \delta_{v, \epsilon}= \epsilon \| (u(x) + \delta_u(x))/(v(x)(v(x) +\delta_v(x)))\bmod x^{m+n+1}\| \end{equation}(for any fixed norm). To estimate the resulting perturbation of the matrix $H_n$, we recall the defintion and some properties of the {\it operator matrix norms} $\|. \|_h$, which are also called {\it subordinate} to and {\it consistent} with the vector norms $\|. \|_h$ of (\ref{equ2.1}) (cf. \cite{gl96}, pp.55--57).  For a matrix $A=(a_{i,j})$, we have $$ \| A \|_h = \sup_{\vec{v}\neq \vec{0}} \|A\vec{v}\|_h/ \|\vec{v}\|,\  h=1, 2, \infty, $$\begin{equation}\label{equ4.3}  \| A \|_1 = \| A^T \|_{\infty} =\max_j \sum_i |a_{i,j}|, \end{equation} $$ \| A \|_2^2 \leq \| A \|_1  \| A \|_{\infty} ,$$ so that \begin{equation}\label{equ4.4} \| A \|_2 \leq \| A \|_1 = \| A \|_{\infty} \end{equation} for a symmetric matrix $A$, in particular for $A=H_n$. (\ref{equ4.3}) implies that the perturbation $\Delta_{\epsilon}$ of $A$caused by the $\epsilon$-perturbations of $u(x)$ and $v(x)$ satisfies $$  \| \Delta_{\epsilon}\|_1 =\| \Delta_{\epsilon}\|_{\infty} \leq \delta_{u, \epsilon} + \delta_{v, \epsilon}$$ for $\delta_{u, \epsilon}$ and $ \delta_{v, \epsilon}$ of  (\ref{equ4.1}), (\ref{equ4.2}). (The equation in the above follows from (\ref{equ4.3})since $\Delta_{\epsilon}$ is a Hankel and therefore symmetric matrix.) Due to (\ref{equ4.4}), we also obtain that \begin{equation}\label{equ4.5} \|\Delta_{\epsilon}\|_2 \leq \delta_{u, \epsilon} + \delta_{v, \epsilon}. \end{equation} To control the resulting impact  of the$\epsilon$-perturbations of$u(x)$ and$v(x)$ on the $\epsilon$-gcds and in particular on their degree $d_{\epsilon}$, onemay  apply some known results on the eigendecomposition of a real symmetric matrix (see \cite{gl96}, Theorems 8.1.1,  8.1.13).\begin{proposition}\label{pro4.5} \begin{description} \item[a)] An $n \times n$ real symmetric matrix $A$ haseigendecomposition \begin{center}$   A=U \bigwedge V^T , $\end{center} where $U$ and $V$ are $n \times n$ orthonormal (unitary) matrices such that %%%%\begin{equation}\label{equ4.4} $$U^TU =V^TV= I, $$%%\end{equation} $\bigwedge = {\rm \  diag}(\lambda_1, \ldots, \lambda_n)$ is a diagonalmatrix, $\lambda_1, \ldots, \lambda_n$ are the eigenvalues of $A$, $\lambda_1 \geq \ldots \geq  \lambda_n$. \item[b)] rank $A = n-z$, where $z$ is the number of the zero eigenvaluesof $A$. \item[c)] For any real symmetric matrix $\Delta$,rank $(A+\Delta) \geq n-z_{\Delta}$, where $z_{\Delta}$ is the number of the eigenvalues of $A$ in the interval from$-\| \Delta \|_2 \sqrt{2}$ to $\| \Delta \|_2 \sqrt{2}$.\end{description}\end{proposition}  The next result follows from the bound (\ref{equ4.5}) and  propositions\ref{pro4.3} and \ref{pro4.5} c). \begin{corollary}\label{cor4.1}Let $d_{\epsilon}$ denote the degree of the $\epsilon$-gcds of $u(x)$ and$v(x)$ and let $z_{\epsilon}$ be the number of the eigenvalues of $H_n$in the interval from $-(\delta_{u,\epsilon} +\delta_{v,\epsilon})\sqrt{2}$ to $(\delta_{u,\epsilon} +\delta_{v,\epsilon})\sqrt{2}$, for $\delta_{u,\epsilon}$ and $\delta_{v,\epsilon}$ of (\ref{equ4.1}), (\ref{equ4.2}). Then $d_{\epsilon} \geq n-z_{\epsilon}$. \end{corollary} Given $\delta_{u,\epsilon} + \delta_{v,\epsilon}$, we compute$z_{\epsilon}$ as follows. We first tridiagonalize $H_n$ by Lanczosalgorithm (\cite{gl96}, ch. 9; \cite{bp94}, sect. 2.3, or \cite{par80}),whose computational cost is dominated by the cost of $n-1$ multiplications of$H_n$ by vectors,  bounded  by $O(n^2 \log n)$ ops (seeproposition \ref{pro4.1}). Then all the eigenvalues of $H_n$can be  approximated within an error bound $2^{-b}$ in$O((n \log^2 n) (\log^2 n +\log b)$ ops \cite{bp91}. Due to corollary \ref{cor4.1}, however, $z_{\epsilon}$can be computed even faster, in $O(n)$ ops, by using the two Sturmsequences formed by the values at $x= -(\delta_{u,\epsilon} +\delta_{v,\epsilon})\sqrt{2}$ and $x = -(\delta_{u,\epsilon} +\delta_{v,\epsilon})\sqrt{2}$ of the characteristic polynomials$p_k(x)$ of the matrices $H_k$ for $k=1, \ldots,n$ (\cite{bp91}, appendixC). \begin{remark}\label{rem4.2} In \cite{cgtw95} and \cite{egl97}, a similar approach wasproposed based on using the customary SVD decomposition of theSylvester and subresultant matrices (cf. our appendix \ref{appa.3}). Since Hankel matrices are symmetric, it is sufficient to usetheir eigendecompositions, thus avoiding the extrasymmetrization step involved in computing the SVD. Decreasing the computational cost by exploiting matrix structure and by only counting thesingular values or eigenvalues (instead of their explicitcomputation) is also a new feature of our algorithm versus\cite{cgtw95} and \cite{egl97}. Such an approach turnedout to be highly effective also for solving a system ofmultivariate polynomial equations \cite{ep97}.  \end{remark} \section{Root neighborhoods and bounding the degree of$\epsilon$-gcds.}\label{sec5} \setcounter{equation}{0} The useful information provided by corollary \ref{cor4.1}has a limitation: the corollary does not promise any non-trivial upper   bound on $d_{\epsilon}$. Indeed, larger eigenvalues of$H_n$ cannot vanish because of the $\epsilon$-perturbations of  $u(x)$ and $v(x)$, but generally we  cannot foresee how many non-zero eigenvalues such a perturbation may turn into 0. In particular, even if we zero one of the absolutely smallestnon-zero eigenvalues in the diagonal matrix of the eigendecomposition of $H_n$ (see proposition \ref{pro4.5}), then  generally such a perturbation would destroy the Hankel structure of $H_n$, andthe resulting matrix would not be associated with any pair ofperturbed polynomials  $u(x)$ and $v(x)$. To remove this limitation, characterizing also  other known approaches to the computation of $\epsilon$-gcds, we will recall the fundamental concept of polynomial root neighborhood, due to \cite{m86}. \begin{definition}\label{def5.1} For a polynomial \begin{equation}\label{equ5.1} p(x) = \sum_{i=0}^k p_i x^i = p_k \prod_{j=1}^k (x-z_j), \  p_k \neq 0, \end{equation} a positive $\epsilon$, and a fixed norm $\|.\|$, we define the$\epsilon$-{\it neighborhood} of $p(x)$, \begin{equation}\label{equ5.2} {\cal N}_{\epsilon} (x) = \{ \tilde{p}(x): \|\tilde{p}(x) - p(x)\| \leq\epsilon, \   \deg \tilde{p}(x) \leq n\},\end{equation} at the $\epsilon$-root-neighborhood of $p(x)$, \begin{equation}\label{equ5.3} Z_{\epsilon} (p) = \{ z, \tilde{p}(z)=0 \  {\rm for \  some} \  \tilde{p}(x) \     {\rm in} \  {\cal N}_{\epsilon} (p) \}.\end{equation} \end{definition} \begin{theorem}\label{the5.1} \cite{m86} Let $Z_{\epsilon} (p)= \bigcup_iC_{\epsilon , i} (p)$, where $C_{\epsilon , i} (p)$ are disjointconnected components of $Z_{\epsilon} (p)$ on the complex plane. Thenevery polynomial $\tilde{p}(x) \in {\cal N}_{\epsilon}(p)$ has thesame number $n_{\epsilon , i} (p)$ of its zeros in $C_{\epsilon , i}(p)$ (counting the zeros with their multiplicity). \end{theorem}Suppose that we have computed a set of disjoint complex domains$D_{\epsilon , i} (u) \supseteq C_{\epsilon , i} (u)$ covering thesecomponents and a similar set of domains $D_{\epsilon , i} (v)\supseteqC_{\epsilon , i} (v) $ for all $i$ and $j$. Let usdefine a bipartite graph $G_{\epsilon}(u,v)$ having two sets of vertices$U =\{u_{i, g}\}$ and $V =\{v_{j, h}\}$, where $n_{\epsilon , i} (u)$ vertices, $u_{i,1},\ldots, v_{i,n_{\epsilon , i} (u)}$, represent$D_{\epsilon , i} (u)$ and $n_{\epsilon , j} (v)$ vertices, $v_{j,1},\ldots, v_{j,n_{\epsilon , i} (v)}$, represent $D_{\epsilon , i} (v)$.  We connect a pair of vertices $u_{i,g}$ and $v_{j,h}$ by an edge if andonly if the domains  $D_{\epsilon , i} (u)$ and $D_{\epsilon , i} (v)$have non-empty intersection. Now we may immediately bound from above the degree $d_{\epsilon}$ of the $\epsilon$-gcds of $u(x)$ and $v(x)$.  \begin{proposition}\label{pro5.1} The degree $d_{\epsilon}$ of an $\epsilon$-gcds of $u(x)$ and $v(x)$cannot exceed the cardinality $d_{\epsilon}^+=|M|$ of a maximum matching$M$  in the associated bipartite graph $G_{\epsilon}(u,v)$. \end{proposition} \section{The complexity of computing an upper bound $d_{\epsilon}^*$ on the degree of$\epsilon$-gcds.}\label{sec6} \setcounter{equation}{0} We recall that the maximum matching in a graph with $n$ vertices can becomputed by using $O(n^{2.5})$ comparisons \cite{hk73}. The domains $D_{\epsilon , i} (p)$ for a given polynomial $p(x)$ of adegree $k$ can be computed in two steps: \begin{description} \item[a)] approximate all the zeros $z_j$ of $p(x)$ within some fixederror bounds $2^{-b_j}$, $b_j=b_j(p, \epsilon)$, \item[b)] apply the known estimates for the upper bounds on theperturbation of $z_j$ caused by the perturbation of $p(x)$ within ${\cal N}_{\epsilon} (p)$. \end{description} Assuming standard and simple normalization of$p(x)$, which brings all its zeros inside a unit disc, the knownalgorithms enable us to perform stage a) by using order of $k$ ops, up toa polylogarithmic factor in $k$ and $b$; furthermore, fully efficient NCparallelization of these algorithms is available \cite{p95}, \cite{p96},\cite{p97}.  On the other hand, the perturbation of azero $z_j$ of $p(x)$ is bounded from above by \begin{equation}\label{equ6.1} \delta_j =(h ! \max_{\tilde{p} \in {\cal N}_{\epsilon}(p)}|\tilde{p}(z_j)/ p^{(h)}(z_j)|+ O(\epsilon^2))^{1/h} \end{equation} provided that $p^{(i)}(z_j)=0, \  i=0, 1, \ldots, h-1$,  and $p^{(h)}(z_j) \neq 0$ (cf. \cite{a80}, p. 82). We substitute into (\ref{equ6.1}) thebound $ |\tilde{p}(z_j)| \leq \epsilon \sum_{i=0}^k |z_j|^i=(1-|z_j|^{k+1})\epsilon / (1-|z_j|)$ and compute someupper bounds $\delta_j^+ $ on$\delta_j$ for all $j$ in $O(k^2)$ ops. (If all zeros of $p(x)$ are simple, we may compute$\delta_j^+$ for all $j$ by using $O(k \log^2 k)$ ops.) The presented approach and the complexity estimates can be extended to the case where  we define the $\epsilon$-neighborhoods $N_{\epsilon, rel}$, and consequently, also the $\epsilon$-root-neighborhoods $Z_{\epsilon, rel}$ based onrelative, rather than absolute error bounds: $$ N_{\epsilon, rel}= \{ \tilde{p}(x)=\sum_{i=0}^k \tilde{p}_i x^i, \  \  |\tilde{p}_i/p_i -1| \leq \epsilon, \  i=0,1,\ldots, k \} $$ for$p(x)$ of (\ref{equ5.1}). (Note that $\tilde{p}_i=0$ if $p_i=0$ under thisdefinition.) The following nontrivial and {\it explicit} (rather than asymptotic in$\epsilon$) {\it upper bounds} on the maximum absolute  perturbation ofthe zeros$z_j$ in  the $N_{\epsilon, rel}$ are due to D. Bini (cf. \cite{bpa}): \begin{proposition}\label{pro6.1}\begin{description} \item[a)] Let $\tilde{p} (x) \in N_{\epsilon, rel}(p)$ for $p(x)$ of (\ref{equ5.1}).Let $z_j \neq 0$ be a zero of $p(x)$ of multiplicity $h=h_j$. Then thereexists $h_j$ zeros$\tilde{z}_j$ of $\tilde{p}(x)$ (counting them with their multiplicity)satisfying the bound \begin{equation}\label{equ6.2}| 1- \tilde{z}_j/ z_j| \leq \epsilon \left( (h ! \left( \begin{array}{c} k \\ h \end{array} \right) (1+\epsilon)\frac{\sum_{i=0}^k |p_i| |z_j|^i }{|z_j| |p^{(h)} (z_j)|} \right)^{1/h},\  \   j=1, \ldots, k.\end{equation}\item[b)] Given $z_1, \ldots, z_k$, the right-hand sideexpressions above can be computed for all $j$  by using $O(k^2)$ ops (orjust $O(k \log^2 k)$ ops if $p'(z_j) \neq 0$ for all $j$). \end{description} Given some upper bounds on the maximum absolute perturbations of thezeros of $u(x)$ and $v(x)$ within $Z_{\epsilon}(u)$ and $Z_{\epsilon}(v)$,respectively, it suffices to use $O(n^2)$ ops to computethe bipartite graph $G_{\epsilon}(u,v)$.  \end{proposition} \section{Extensions to computing $\epsilon$-gcds.}\label{sec7} \setcounter{equation}{0} Having computed a maximum matching $M$ in the bipartite graph$G_{\epsilon} (u,v)$, we may select $d_{\epsilon}^+= |M|$ points $w_k$ in the $d_{\epsilon}^+$ intersections $D_{\epsilon, i(k)}(u) \bigcap D_{\epsilon, j(k)}(v), \   k=1, \ldots,d_{\epsilon}^+$, of the pairs of the domains associated with$d_{\epsilon}^+ = |M|$ edges of $M$. The choice of the points $w_k$ in thedomains is not unique, and various heuristics could be proposed; forinstance, we may choose $w_k=(u_k +v_k)/2$, where $u_k$ and $v_k$ denotethe computed approximations to the zeros of $u(x)$ and $v(x)$ lying inthe domains $D_{\epsilon, i(k)}$ and  $D_{\epsilon, j(k)}$,respectively. Once the values $w_1, \ldots, w_{d_{\epsilon}^+}$ are available, wecompute the candidate $\epsilon$-gcd, $\prod_{k=1}^{d_{\epsilon}^+}(x-w_k)$, and test if it a common $\epsilon$-divisor of $u(x)$ and$v(x)$. If so, it is output as an $\epsilon$-gcd. Otherwise we may trydistinct heuristics for defining $\epsilon$-neighborhoods and/orchoosing $w_k$, for $k=1, \ldots, d_{\epsilon}$. If this does not worktoo, we may increase $\epsilon$ or output FAILURE and try anotherapproach, where we may use the upper bound $d_{\epsilon}^+$. \section{Approximating the gcds of polynomials represented by theirzeros.}\label{sec8} \setcounter{equation}{0} In the beginning of this paper, we have shifted from gcds to$\epsilon$-gcds to escape from the unstable dependence of the degree ofthe gcd on the perturbation of the input coefficients. The next simpleexample shows that the remedy does not work perfectly well since thedegree of $\epsilon$-gcds is also very sensitive to the smallperturbations of these coefficients, at least in the case where the inputpolynomials have multiple or clustered zeros. {\bf Example 8.1.} Let$n$ be even and large,$\epsilon =2^{-n}$, $u(x)=x^n$,$v(x)=(x-1)^n$. Then gcd $(u(x),v(x))=1$. Moreover,the distance between the only zero$x=0$ of $u(x)$ and the only zero $x=1$ of $v(x)$ is large enoughto suggest that, under any reasonable definition, an approximate gcdof $u(x)$ and $v(x)$ should be a constant too. $\epsilon$-gcds of $u(x)$ and $v(x)$, however, have degree 1: forinstance, we have an $\epsilon$-gcd  $x-0.5$, which divides bothpolynomials $u^* (x)=u(x)-2^{-n}$ and $v^* (x)=v(x)- 2^{-n}$ satisfying (\ref{equ1.1}) and (\ref{equ1.2}) for $\epsilon=2^{-n}$.  As a way out of this trap, we propose to  define the $\epsilon$-neighborhoods based on the perturbation of thepolynomial zeros, rather than the coefficients. We will study such a natural direction in this section and in the next one and will showcomputational efficiency of this approach. To distinguish this study fromthe previous definitions, we will talk about$\delta$-gcds, rather than $\epsilon$-gcds, and will rely on therepresentation of the input polynomials $u(x)$ and$v(x)$ of (\ref{equ2.2}) by their zeros: \begin{equation}\label{equ8.1} u(x)= u\prod_{i=0}^m (x-y_i), \   v(x)= v\prod_{i=0}^m (x-z_i), \   \      uv\neq 0. \end{equation}The $\delta$-neighborhoods of $u(x)$ and$v(x)$ will be defined as the sets of polynomials $N^*_{\delta}(u)=\{ u^*(x)= u\prod_{i=0}^m (x-y_i^*) \}$ and  $N^*_{\delta}(v)=\{ v^*(x)= v\prod_{j=0}^m (x-z_j^*) \}$ with $y_i$ and $z_j$ satisfying $$  |y_i^* - y_i| \leq \delta,   \    |z_j^* - z_j| \leq \delta \   {\rmfor \  all} \  i \  {\rm and} \  j.$$  The common monic divisors of all pairs $u^*(x)$,$v^*(x)$ from these  $\delta$-neighborhoods will be called $\delta$-{\itdivisors} of $u(x)$ and $v(x)$, and such $\delta$-divisors of the maximumdegree $d_{\delta}$ will be called $\delta$-{\it gcds} of $u(x)$ and$v(x)$.  To compute a $\delta$-gcd, we proceed as in sections \ref{sec5}--\ref{sec7}. The bipartite graph $G$ is defined by its two vertex sets, $Y=\{y_1, \ldots, y_n \}$ and $Z=\{z_1, \ldots, z_n \}$ with pairs $y_i$ and $z_j$ of vertices connected by an edge if only if $ |y_i - z_j| \leq 2 \delta$. A maximum matching $M = \{ (y_{i_1}, z_{j_1}), \ldots, (y_{i_r}, z_{j_r})\}$ in $G$ immediately defines a correct $\delta$-gcd  $g^*(x)$ of $u(x)$ and $v(x)$, $g^*(x)=\prod_{q=1}^r(x - (y_{i_q}+z_{j_q})/2)$. The computation involves $O(n^2)$ ops and comparisons to construct thegraph $G$ and $O(n^{2.5})$ comparisons to compute a maximum matching in$G$. It takes extra $O(n \log^2 n)$ ops to compute thecoefficients of $g^*(x)$ (\cite{bp94}, pp. 25). The next proposition shows that the polynomial $g^*(x)$ is also acommon $\epsilon$-divisor of  $u(x)$ and $v(x)$ if $\delta \leq(\epsilon)^{1/n} -1$. \begin{proposition}\label{pro8.1} Let \begin{equation}\label{equ8.2} \delta\leq  (1 + \epsilon)^{1/n}-1 \ . \end{equation} Then any $\delta$-gcd of the polynomials $u(x) $ and $v(x)$ of (\ref{equ8.1}) is the gcd oftwo polynomials$u^* (x)$ and $v^* (x)$ satisfying  (\ref{equ1.1}) and(\ref{equ1.2}),  provided that $m \leq n$.\end{proposition}{\bf Proof.} The value$\parallel\tilde{u}(x)-u(x)\parallel/ \parallel u(x)\parallel $ reaches its maximum for \\$\displaystyle{ u(x)=u\prod_{i=1}^m (x-y_i ), \tilde{u}(x) = u\prod_{i=1}^m (x-y_i -\delta ) }=u(x-\delta)$and for some nonnegative $y_i$,$i=1,\ldots ,n$. For such $u(x)$ and $\tilde{u}(x)$, we have$$\parallel \tilde{u}(x)-u(x)\parallel =\parallel\sum_{i=1}^mu^{(i)} (x) \delta^i /i!\parallel\leq\sum_{i=1}^m \parallelu^{(i)} (x)/i!\parallel\delta^i \leq$$$$\parallel u(x)\parallel\sum_{i=1}^m \delta^i \left ( \begin{array}{c}m \\ i \end{array} \right )= \parallel u(x)\parallel((1+\delta)^m -1 ) \ .$$If $m \leq n$ and (\ref{equ8.2}) holds,then $(1+\delta)^m \leq (1+\delta)^n \leq 1+ \epsilon $, and consequently,$$\parallel \tilde{u}(x)-u(x)\parallel \leq 2^{-b} \parallel u\parallel \ .$$Similarly, we deduce that $\parallel\tilde{v}(x)-v(x)\parallel \leq \epsilon \|v(x)\|$ if $|z_j-\tilde{z}_j| \leq \delta$, $j=1, \ldots, n$, for $\delta$ satisfying (\ref{equ8.2}).  \hfill $\Box$\section{Simplified computation of $\delta$-gcds.}\label{sec9} \setcounter{equation}{0} Let us simplify the computation of a maximum matching by using {\itelastic} $\delta$, that is, by allowing to adjust $\delta$ by increasingit dynamically. Our first observation is that a maximum matching in $G$ isimmediately available if we know that all the componente of $G$ arebipartite cliques. Now, there are at most $c \leq 2m-1$ edges in any simple path in $G$, andeach edge connects two points $y_i$ and $z_i$, where $|y_i-z_i| \leq2\delta$. If we replace $\delta$ by $\delta_1=(4n -2)\delta$, then it can be immediately estimated that everycomponent of $G$ will become a bipartite clique in the new bipartitegraph $G_1$. On the other hand, if the minimum bipartite distance betweenthe pair of the components of $G$ exceeds $\delta_1$, then none of thepairs of the components of $G$ can be connected by an edge in $G_1$, and therefore, all components of $G_1$ are bipartite cliques; inthis case  we easily compute a maximum matching in $G_1$ and a$\delta_1$-gcd. Otherwise, $G_1$ has fewer components than $G$. In thiscase we may apply the same techniques recursively, replace $\delta$ by$\delta_q=(4n -2)^q\delta$ and compute bipartite graphs $G_q$ for $q=1,\ldots, I$. We will stop where the graph $G_I$ has as many components as$G_{I-1}$; then all of these components must be bipartite cliques, and wewill immediately output a $\delta_I$-gcd. Clearly, $I \leq m$, andactually $I \leq m-1$ since we may ignore thecomponents-singletons. The entire computation reduces to computing the connected components in$G=G_0$ (which requires $O(mn)$ comparisons \cite{clr90}) and updatingthem in $G_q$ for $q=1, \ldots, I$. For updating, one may create and recursively update the list of the minimum distances from each vertex $y_i$ to thevertices $z_j$ in each component of $G_q$. The entire recursive updatingcan be performed in $O(mn)$ comparisons and link operations, so thatwe have the bound $O(mn)$ on the entire cost of computing a$\hat{\delta}$-gcd, for  $\hat{\delta}=\delta_I$, $I \leq m-1$, $m \leq n$. This algorithm can beimplemented in $O(\log^2n)$ parallel time by using $O(mn)$ work. To achieve this, we apply binary search in the interval $[0, m-1]$ tofind $I$ and at each step use the known effective parallel algorithms for  connected components (cf. \cite{hz96}). Note that parallel computation of a maximum matchingis a harder  problem \cite{gp88}.\section{Extension to $\epsilon$-Pad\'{e} approximation and $\delta$-Pad\'{e} approximation.}\label{sec10} \setcounter{equation}{0} Our study of $\epsilon$-gcds and $\delta$-gcds can be extended to thestudy of $(k, l)$-th Pad\'{e} approximations allowing inexact input, thatis, inexact coefficients of the formal power series $h(x)$.  Byfollowing  the patterns of section~\ref{sec3}, we may compute the $(k-d, l-d)$-th Pad\'{e} approximations $s_d(x)$, $t_d(x)$ of $h(x)$ for $d=0,1, \ldots, \min(k,l)$, then test if $$ \| (h(x) - s_d(x)/t_d(x)) \bmod x^{k+l-2d} \| \leq \epsilon$$ for a fixed $\epsilon$, and output the pair $s_d(x)$, $t_d(x)$ that passed the test and has the greatest degree $d$. Alternatively, we maycompute only the single $(k, l)$-th Pad\'{e} approximation  $s(x)$, $t(x)$ of $h(x)$ and then apply the techniques of sections\ref{sec3}, \ref{sec5}--\ref{sec9} to compute$g(x)$, an $\epsilon$-gcd (or $\delta$-gcd) of $s(x)$ and $t(x)$. Finally, we may compute $\epsilon$-quotients (or $\delta$-quotients) for $s_d(x)=s(x)/g(x)$ and  $t_d(x)=t(x)/g(x)$ and test if the norm $\|s_d(x)- h(x)t_d(x)\|$ is sufficiently small.\section{Discussion.}\label{sec11} \setcounter{equation}{0} Choosing an appropriate measure for the distances in the study ofapproximate gcds is a critical but still largely open problem. Usingrelative errors of the perturbation and/or transition to measuring theperturbation in terms of the zeros of $u(x)$ and $v(x)$ seem to beappropriate  directions, partly validated in this paper. Further progress in our understanding and computing $\epsilon$-gcdsdepends also on further analysis of $\epsilon$-root-neighborhoods anddevising effective techniques for their computation. Some improvement of the computation of $\delta$-gcds may rely on thesimplification of approximating the zeros of $u(x)$ and $v(x)$: since weonly care about nearly common zeros, we may seek effective methods fordiscarding the factors of $u(x)$ and $v(x)$ whose zero sets do not comeclose to each other. Then the deflation of $u(x)$ and $v(x)$ woulddecrease the input size. \begin{thebibliography}{92} \bibitem[A80]{a80}K. E. Atkinson, {\it Introduction to Numerical Analysis}, Wiley and Sons, New York, 1980.\bibitem[B81]{b81}  S. Barnett, Generalized Polynomials and LinearSystems  Theory, {\it Proc. 3rd IMA Conf. Contr. Theory Edit.}, 3--30, 1981.\bibitem[BGY80]{bgy80}R.B. Brent, F.G. Gustavson, D.Y.Y. Yun, Fast Solution of Toeplitz Systemsof Equations and Computation of Pad\'{e} Approximations, {\it J.of Algorithms}, 1, 259--295, 1980.\bibitem[BP91]{bp91}   D. Bini, V. Y. Pan,  Parallel Complexity of TridiagonalSymmetric  Eigenvalue Problem, {\it Proc. 2nd Ann. ACM-SIAM Symposiumon Discrete Algorithms,} 384--393, ACM Press, New York, and SIAMPublications,  Philadelphia, Pennsylvania, 1991.\bibitem[BP94]{bp94}D. Bini, V.Y. Pan, {\it Polynomial and Matrix Computations,vol. 1: Fundamental Algorithms}, Birkhaeuser, Boston, Massachusetts,1994. \bibitem[BP,a]{bpa}D. Bini, V.Y. Pan, {\it Polynomial and Matrix Computations,vol. 2}, Birkhaeuser, Boston, Massachusetts,to appear.\bibitem[CLR90]{clr90}T. H. Cormen, C. E. Leiserson, R. L. Rivest, {\it Introduction to Algorithms}, MIT Press, Cambridge, Mass., and McGraw-Hill, New York, 1990. %%\bibitem[CW90]{cw90} D. Coppersmith, %%S. Winograd, Matrix Multiplication via %%Arithmetic Progressions, {\it J.%%of Symbolic Computation}, 9, 3, 251-280, 1990.\bibitem[CGTW95]{cgtw95}R.M. Corless, P.M. Gianni, B.M. Trager, S.M. Watt,The Singular Value Decomposition for Polynomial Systems, {\it Proc. Intern.Symp. on Symbolic and Algebraic Comp. (ISSAC '95)}, 195-207, ACMPress, New York, 1995.%%%\bibitem[D83]{d83}  J. D. Dixon, %%Estimating Extremal Eigenvalues and%%Condition Numbers of Matrices, {\it SIAM J. %%of Numerical Analysis}, 20,4, 812--814, 1983.\bibitem[EGL97]{egl97}I. Z. Emiris, A. Galligo, H. Lombardi, Certified ApproximatePolynomial GCDs, to appear in {\it J. of Pure and Applied Algebra}, 1997.\bibitem[EP97]{ep97} I. Z. Emiris, V. Y. Pan, The Structure of Sparse Resultant Matrices, {\it Proc. of Annual ACM-SIGSAM International Symp. on Symb. and Alg. Comp. (ISSAC'97)}, ACM Press, New York, 1997. %%\bibitem[G84]{g84}J. von zur Gathen, %%Parallel Algorithms for Algebraic 	Problems, %%{\it SIAM J. on Computing}, 13, 4, 802-824, 1984.\bibitem [G95]{g95}A. Galligo, On Some Difficulties of Approximating  Polynomial GCD, {\it Lecture at AMS-SIAM Conference: Mathematics of Numerical Analysis: Real Number Algorithms,} Park City, Utah, July-August 1995. \bibitem [GL96]{gl96}   G.H. Golub, C.F. Van Loan, {\it Matrix Computations,} Johns  Hopkins Univ. Press,  Baltimore,  Maryland,  1996.\bibitem[GP88]{gp88}Z. Galil, V. Y. Pan, Improved Processor Bounds for CombinatorialProblems in RNC, {\it Combinatorica}, 8, 2, 189-200, 1988.\bibitem[Gr72]{gr72}W.B. Gragg, The Pad\'{e} Table and Its Relation to Certain Algorithmsof Numerical Analysis, {\it SIAM Review}, 14, 1, 1-62, 1972.\bibitem[HK73]{hk73}J.E. Hopcroft, R.M. Karp, An $n^{5/2}$ Algorithm for MaximumMatching in Bipartite Graphs, {\it SIAM J. on Computing}, 2,225-231, 1973.\bibitem[HS95]{hs95}V. Hribernig, H.J. Stetter, Detection and Validation ofClusters of Polynomial Zeros, {\it J. Symb. Comp.}, 1995.\bibitem [HZ96]{hz96} S. Halperin, U. Zwick, Optimal Randomized EREW PRAM Algorithms for Finding Spanning Forests and for Other Basic Graph Connectivity Problems, {\it Proceedings of the 7th Ann. ACM-SIAM Symposium on DiscreteAlgorithms},  384--393, ACM Press, New York, and SIAM Publications, Philadelphia, Pennsylvania, 1996.%%\bibitem[II86]{iis86}A. Israeli, A. Itai, %%A Fast and Simple Randomized Parallel %%Algorithm for Maximal Matching, {\it %%Information Proc. Letters}, 22, 77-80, 1986.%\bibitem[IS86]{is86} A. Israeli, %Y. Shiloach, An Improved Parallel Algorithm %for Maximal Matching, {\it Information Proc. %Letters}, 22, 57-60, 1986.\bibitem[J92]{j92}J. J\'{a}J\'{a}, {\it An Introduction to Parallel Algorithms},Addison-Wesley, Reading, Mass., 1992.\bibitem[K66]{k66} W. Kahan, Numerical Linear Algebra, {\it Canadian Math. Bull.}, 9,757--801, 1966.\bibitem[K80]{k80} T. Kailath, {\it Linear Systems}, Prentice-Hall, Englewood Cliffs, New Jersey, 1980.\bibitem[K87]{k87} N. Karcanias, Invariance Properties and Characterisations of the Greatest Common Divisor of a Set of Polynomials, {\it Int. Journal of Control}, 40, 1751--1760, 1987. \bibitem[KL96]{kl96}N. Karmarkar, Y. N. Lakshman, Approximate Polynomial Greatest CommonDivisor and Nearest Singular Polynomials, {\it Proc. Ann. ACM-SIGZAMIntern. Symp. on Symb. and Algebraic Comp.,} 35--39, 1996. \bibitem[KM94]{km94}N. Karcanias, M. Mitrouli, A Matrix Pencil Based Numerical Methodfor the Computation of GCD of Polynomials, {\it IEEE Trans. onAutomatic Control}, 39, 5, 977-981, 1994.\bibitem [L81]{l81}  D. Lazard, Resolution des Systems d'EquationsAlgebrique, {\it Theoretical Computer Science, } 15, 77--110, 1981. \bibitem[M86]{m86} R. G. Mosier, Root Neighborhoods of Polynomials, {\itMath. of Computations}, 47, 175, 265--273, 1986. \bibitem[NS91]{ns91}M.-T. Noda, T. Sasaki, Approximate GCD and Its Application toIll-Conditioned Algebraic Equations, {\it J. Computational and AppliedMathematics}, 38, 335-351, 1991.\bibitem[P92]{p92}V.Y. Pan, Complexity of Computations with Matrices andPolynomials, {\it SIAM Review}, 34, 2, 225-262, 1992.\bibitem[P95]{p95}V.Y. Pan,Optimal (up to Polylog Factors) Sequential and Parallel Algorithmsfor Approximating Complex Polynomial Zeros, {\it Proc. 27th Ann. ACMSymposium on Theory of Computing}, 741-750, ACM Press, New York, 1995.\bibitem[P96]{p96}V.Y. Pan,Optimal and Nearly Optimal Algorithms for Approximating PolynomialZeros, {\it Computers and Mathematics (with Applications)}, 31, 12, 97-138, 1996.\bibitem[P96,a]{p96a}V. Y. Pan, Numerical Computation of a Polynomial GCD and Extension, Research Report 2969,{\it INRIA, Sophia-Antipolis, France}, 1996.\bibitem[P97]{p97}V.Y. Pan, Solving a Polynomial Equation: Some History and RecentProgress,  {\it SIAM Review}, 39, 2, 1997. \bibitem[Par80]{par80}B. N. Parlett, {\it The Symmetric Eigenvalue Problem,} Prentice-Hall,Englewood Cliffs, New Jersey, 1980. \bibitem [PB73]{pb73} J. S. Pace, S. Barnett, Comparison of Algorithms for Calculation of GCD of Polynomials, {\it Int. Journal Systems Sci.}, 4, 211--223, 1973. \bibitem[Sc85]{sc85}A. Sch\"{o}nhage, Quasi-GCD Computations, {\it J. of Complexity},1, 118-137, 1985.\end{thebibliography}\appendix \section{Some known approaches to $\epsilon$-gcds.}\label{appa}\setcounter{equation}{0}  \subsection{The quadratic programming approach.}\label{appa.1}   Seeking a linear $\epsilon$-divisor of two polynomials $u(x)$ and $v(x)$ of (\ref{equ2.2}), we may minimize the 2-norms, $$ \| u(x) -(x -\alpha) \Phi (x) \|_2 \  {\rm and } \   \| v(x) -(x -\alpha) \Psi (x) \|_2 ,$$  over all $\alpha$ and polynomials $\Phi (x)$ and $\Psi (x)$ of degrees $m-1$ and $n-1$, respectively, and then check the $\epsilon$-divisibility of $u(x)$ and $v(x)$ by $x-\alpha$. This procedure was proposed and extended to nonlinear $\epsilon$-divisors in \cite{kl96}. The computation is reduced toquadratic programming, and the estimated number of ops involved ispolynomial in$n$  (of unspecified degree) and exponential in $d$. \subsection{The Euclidean algorithm approach. }\label{appa.2}%to gcds and  $\epsilon$-gcds} \setcounter{equation}{0}  To perform Euclidean algorithm, we write $r_0(x)=v(x)$, $r_1(x)=u(x)$ and then recursively compute $r_{i+1} (x)=r_{i-1}(x) \bmod r_i(x)$, $i=1, 2, \ldots, l$, until wearrive at  $r_{l+1} (x)=0$. Then we output the exact gcd, $r_l (x)= g(x) = {\rm gcd} (u(x), v(x))$. The algorithm has a fastimplementation in $O(n \log^2 n)$ ops (\cite{bp94}, pp. 37--39). To compute an $\epsilon$-gcd, one may stop the computations earlier, at the smallest $k$ for which $r_k(x)$ is a common $\epsilon$-divisor %ibility of $u(x)$ and $v(x)$ \cite{sc85}, \cite{egl97}. The algorithm does not ensure that the degree of the computed common $\epsilon$-divisor is maximum (see a counter-example in \cite{egl97}) butproduces a lower bound on this maximum.  \subsection{The subresultant approach.}\label{appa.3} \setcounter{equation}{0}  For two polynomials $u(x)$ and $v(x)$ of (\ref{equ2.2}), their {\itextended gcd} is the triple $(g(x), \beta(x), \gamma(x))$ formed by the gcdand the two cofactors satisfying the following relations: \begin{equation}\label{equa.1}\beta(x) u(x) + \gamma(x) v(x) = g(x), \end{equation}\begin{equation}\label{equa.2} g(x) = {\rm gcd}(u(x), v(x)), \  \  \deg g(x)=d , \end{equation}\begin{equation}\label{equa.3} \deg \beta(x) < n-d, \  \   \deg \gamma(x) < m-d.\end{equation} Clearly, the polynomial $g(x)$ can be recovered from (\ref{equa.1}) if weknow the cofactors $\beta(x)$ and $\gamma(x)$. Their coefficient vectors $\vec{\beta}$ and $\vec{\gamma}$, ofdimensions $n-d$ and $m-d$, respectively, satisfy the vector equation \begin{equation}\label{equa.4}S_k \left[ \begin{array} {c} \vec{\beta}\\  \\  \\ \vec{\gamma} \end{array}     \right]= \left[ \begin{array} {c} 0\\ \vdots  \\ 0 \\  1 \end{array} \right], \end{equation}where $k=d$ and $S_k=S_k(\vec{u},\vec{v})$ is called a {\it subresultant  matrix}. \begin{proposition}\label{proa.1}The system (\ref{equa.4}) isnon-singular for$k=d$, and is singular for $k < d$. \end{proposition} If $d$ is known, we may compute the vectors $\vec{\beta}$ and$\vec{\gamma}$ by solving the linear  system (\ref{equa.4})and then compute$g(x)$, the gcd, from (\ref{equa.1}). The same algorithm can be applied to compute a candidate $\epsilon$-gcd,$g_{k, \epsilon}(x)$, for a  fixed degree $k=d_{\epsilon}$; atthe end we may test the $\epsilon$-divisibility of $u(x)$ and $v(x)$ by thecomputed polynomial $g_{k, \epsilon}(x)$ (see proposition~\ref{pro2.1}). Even if the degree $d=d_{\epsilon}$ of $\epsilon$-gcd is notavailable, we may make a search for the maximum $k$, $k \leq m$, forwhich the system (\ref{equa.4}) has a solution such that the associatedpolynomial $g(x)$ is a common $\epsilon$-divisor of $u(x)$ and$v(x)$. Such a search can be narrowed by discarding those values $k$ forwhich the $(m+n-k+1)$-st largest singular value of $S_0$ exceeds $\epsilon$ (see \cite{cgtw95}, \cite{egl97}). The matrices $S_k$ havestructure that is similar to one of Hankel matrices $H_k$ but a littlemore complicated, in particular the matrices $S_k$ are not symmetric. Onthe other hand, it is easier to estimate the perturbations of $S_k$(rather than ones of $H_k$) caused by the perturbations of $u(x)$ and$v(x)$.\section{Verification of the $\epsilon$-divisibility under aweighted 2-norm.}\label{appb}\setcounter{equation}{0}  Let us prove proposition \ref{pro2.2}. The equation $p(x)-d(x)g(x) =0$ isequivalent to the linear system $$ H(\vec{d}) \vec{g}_r= \vec{p}, $$where $\vec{g}_r=(g_{k-l}, \ldots, g_0)^T$ is the reverse of the vector $$\vec{g}=(g_0, \ldots, g_{k-l})^T$$and $H(\vec{d})$ is a $(k+1) \times (k+1-l)$ Hankel matrix of full rank(cf. \cite{bp94} or \cite{cgtw95}. Minimizing $\| p(x) - d(x) g(x)\|_{2, \vec{w}}$ is equivalent to finding aweighted least-squares solution to the latter linear system, which isthe solution to the non-singular linear system of $k-l +1$ normalequations \cite{a80}, $$ H^T(\vec{d}) \  {\rm  diag} \   (\vec{w}) H(\vec{d}) =H^T(\vec{d})\vec{p}, $$ where diag $(\vec{w}) $ is the  $(k-l +1) \times (k-l+1)$ diagonal matrixwith the diagonal entries $w_1, \ldots, w_{k-l+1}$. This system can be solved by the conjugate gradient method,\cite{gl96}, section 10.2, at the cost of performing $O((k-l) \mu)$ ops, $\mu$ denoting the number of ops formultiplication of the matrix $H^T(\vec{d})$ diag$(\vec{w})H(\vec{d})$ by a vector, $\mu = O(k \log (k-l))$, byproposition \ref{pro4.1}. \hfill $\Box$ \end{document} \section{Matrix norms and singular values.}\label{sec4}\setcounter{equation}{0}  By using the $h$-norms of vectors, we define the subordinate (operator) $h$-norms of matrices, $$ \|A\|_h = \sup_{\vec{v} \neq \vec{0}} \|A \vec{v}\|_h/\|\vec{v}\|_h.$$In some applications to the $\epsilon$-gcds, we will need to estimate the distance from a fixed matrix $A$ to the closest singular matrix, which can be simply expressed via $\|A^{-1}\|_h$ if $A$ is nonsingular (cf. \cite{k66}): \begin{equation}\label{equ4.1} d_h(A)=\min_{\det B =0} \|A-B\|_h = 1/\| A^{-1}\|_h.\end{equation}As soon as $A^{-1}$ is available, we may compute or estimate $d_h(A)$ for $h=1,2, \infty$ based on the next proposition. \begin{proposition}\label{pro4.1} \cite{gl96}, pp. 56--57. For any $k \times k$ matrix $A=(a_{i,j})$, $i,j=1, \ldots, k$, we have $$1/k \|A\|_h \leq \max_{i,j} |a_{i,j}| \leq  \|A\|_h, {\rm \  for \   } h=1,2,\infty,$$$$ \|A\|_{\infty} = \|A^T\|_1 = \max_{1 \leq i \leq k} \sum_{j=1}^k |a_{i,j}|, $$ $$ \sqrt{k} \|A\|_h \leq  \|A\|_2 = \|A^T\|_2 \leq \sqrt{k} \|A\|_k {\rm \  for} \  h=1 {\rm \  and} \  h=\infty,$$ $$  \|A\|_2^2 \leq  \|A\|_1 = \|A^T\|_{\infty}.$$ \end{proposition} In our study of $\epsilon$-gcds, we will also seek the distances from a fixed matrx $A$ to the closest matrices of a rank at most $r$, for fixed $r<k$, which leads us to the concept of {\it singular values} of a matrix.\begin{definition}\label{def4.1} diag$_k(a_1, \ldots, a_r)$ denotes the $ k \times k$ (diagonal) matrix havingzero entries every where except for the first $r$ entries $a_1, \ldots, a_r$ of its diagonal. $I_k$ is the $ k \times k$ identity matrix, $I={\rm diag}_k(1,\ldots,1)$, $r=k$. A $ k \times k$ real matrix $A$ is called orthogonal (alsoorthonormal and unitary) if $A^TA=I_k$. \end{definition} \begin{definition}\label{def4.2} (SVD) Singular value decomposition of a real $ k \times k$ matrix $A$ is a triple of $ k \times k$ matrices $(U, \Sigma, V)$ satisfying the relations $A=U\Sigma V^T$, $U^TU=I_k$, $V^TV=I_k$, $\Sigma={\rm diag}_k(\sigma_1, \ldots,\sigma_r)$,  $\sigma_1 \geq \ldots \geq \sigma_r >0$, $r=$rank $A \leq k$, $\sigma_i=\sigma_i(A)$ is the $i$-th {\it singular value} of $A$, and the $i$-th columns of $U$ and $V$ are the left and right $i$-th {\it singular vectors}, respectively, $i=1, \ldots, r$. \end{definition}SVD has the following properties (\cite{gl96}), pp. 70--73): \begin{proposition}\label{pro4.2}  Any real matrix $A$ has its SVD. \end{proposition} \begin{proposition}\label{pro4.3}  $\min_{ {\rm rank} \  B < h} \|A-B\|_2 = \sigma_h(A)$, $h \leq {\rm rank} \    A$. \end{proposition} By proposition \ref{pro4.3}, $\sigma_h = \sigma_h(A)$ is the 2-distance from $A$ to the closest matrix of a rank less than $h$, and $\sigma_h$ are the quantitiesthat we will use in our study of $\epsilon$-gcds (we will not use the singularvectors). \section{Approximation of the singular values of a matgrix.}\label{sec5} \setcounter{equation}{0}  To approximate the singular values of a $k \times k$ real matrix  $A=U\Sigma V^T$, we will first shift to the real symmetric matrix $S$, by choosing $S=A^TA$ or $S=AA^T$ and noting that \begin{equation}\label{equ5.1} A^TA= V \Sigma^2 V^T, \  \    AA^T =U \Sigma^2 V^T,\end{equation}then compute a real symmetric tridiagonal matrix $$ T= QQ^T S Q, \  \  Q^TQ=I, $$  so that $T=Y\Sigma^2 Y^T$ for some orthogonal matrix $Y$, and finally approximate the singular values $\sigma_1^2, \ldots, \sigma_r^2$ of $T$ and $\sigma_1, \ldots, \sigma_r$ of $A$. Hereafter, we will  refer to real symmetric matrices as {\it rsmatrices} and will refer to rs tridiagonal matrices as {\it rstmatrices}. The customary transition from the matrix $S$ to $T$ is by Lanczos algorithm\cite{gl96} which uses $O(k^2 +kM(S))$ ops, assuming here and hereafter that {\it the matrix $S$ can be multiplied by a vector in $M(S)$ ops}. More precisely this bound holds if the Lanczos process is not interruptedby its degeneration more than $O(1)$  times. Each interruption is highlyunlikely for any input matrix $S$ due to involving random auxiliary vectorsin the computations. Even in the extremely unlikely) casae of the maximum number of $k-1$ interruptions, $O(k^3$ ops still suffice. Now, given a $k \times k$  rst matrix, we have the following estimates for the costof approximating its singular values: \begin{proposition}\label{pro5.1}\cite{bp94}. The singular values of a $k \times k$  rst matrix $T$ can be approximated within the error bound$\delta \|T\|_1$ at the cost of performing $O((k \log^2k)(\log k + \log^2\log(1/\delta)))$ ops.   \end{proposition} \begin{proposition}\label{pro5.2} Given a $k \times k$  rst matrix $T$ and a positive $\epsilon^2$, it suffices touse $O(k)$ ops to compute the number $k-S(\epsilon^2)$ of singular values $\sigma_k^2, \ldots, \sigma_{S(\epsilon^2)+1}^2$.\end{proposition} {\bf Proof.} Let $p_i(x)$ denote the characteristic polynomial of the $i \times i$ laeding principle (northwestern) submatrix of $T$, $i=1, \ldots, k$. The polynomials $p_1(x), \ldots, p_k(x)$ form a Sturm sequence: the  number of signagreements between the pairs of the sequence of the values $p_i(\alpha), \ldots, p_k(\alpha)$ for any real $\alpha$ equals the number of thesingular values of $T$ that are less than $\alpha$, \cite{par80}, pp. 131; inparticular this holds for $\alpha= \epsilon^2$. The computation of $p_1(\epsilon^2), \ldots, p_k( \epsilon^2)$ takes $O(k)$ ops, due to the simple three-term recurrence expressions of $p_{i+1}(x)$ via $p_i(x)$ and $p_{i-1}(x)$ for all $i$, cf. \cite{bp91}, appendix~\ref{appc}. \hfill $\Box$\section{Randomized approximation of the smallest singular value of a matrix.}\label{sec6}  \setcounter{equation}{0}  Trying to approximate $d_h(A)$ based on (\ref{equ4.1}), one may have practical difficulties with the inversion of $A$ where $A$ is nearly singular matrix. The next approach applies to the matrix \begin{equation}\label{equ6.1} S=A^TA \  \   {\rm or} \  \  S=AA^T \end{equation} (see  (\ref{equ5.1})) and avoids matrix inversion. It relies on the simpleobservation that \begin{equation}\label{equ6.2} \sigma_i( \eta I_k -S) = \eta -\sigma_i(S) \end{equation}  for any $\eta \geq \sigma_1(S)$ and on the following result from \cite{d83}: \begin{proposition}\label{pro6.1} Let $S$ be a $ k \times k$ rs matrix of (\ref{equ6.1}), with the largest singular value $\sigma_1=\sigma_1(S)$. Let $\vec{v}$ be a random vector chosen under the uniform probability distribution on the unit$(k-1)$-dimensional sphere $S_{k-1} = \{ \vec{x}, \| \vec{x} \|=1\}$.Then $\vec{v}^TA^q\vec{v} \leq \sigma_1$, $q=1, 2, \ldots$, and forany positive $\Theta$ any fixed $q$ we have with a probability at most $p_{k, \Theta, q}=0.80 \sqrt{k/(1+\Theta)^q}$, that   $\vec{v}^TA^q\vec{v} > (1+ \Theta) \sigma_1$. \end{proposition} Based on these results, we devise the following algorithm. \begin{algorithm}\label{alg6.1} $$  \   \   $$ {\bf Input:} a $k \times k$ rs matrix $S$ of (\ref{equ6.1}) and two positivescalars, $\mu \leq 1$ and $t$. \\{\bf Output:} $\sigma_k^+$ such that  \begin{description} \item[a)] $\sigma_k(S) - \sigma_k^+ > t$ with a probability at most $\mu$, where $\sigma_k(S)=0$ if rank $S <k$, and \item[b)] $\sigma_k^+ \geq \sigma_k(S)$.\end{description} {\bf Computations.} \begin{enumerate} \itemCompute an upper bound $\sigma_1^+$ on $\sigma_1(S)=\|S\|_2$. \item Compute \begin{equation}\label{equ6.3} q=|\frac{\log (25 k \mu^2/16}{\log (1+t/\sigma_1^+}| .\end{equation}\item Compute a random vector $\vec{v}$ chosen under the uniform probabilitydistribution on the sphere $S_{k-1}$. \item Compute and output the value \begin{equation}\label{equ6.4} \sigma_k^+ = \sigma_1^+ -\vec{v}^T (\sigma_1^+ I -S)^q \vec{v}.\end{equation} \end{enumerate} \end{algorithm} To prove {\it correctness} of the algorithm, first obtain from (\ref{equ6.2}) that $$ \sigma_1^+  - \sigma_k(S) = \sigma_1(\sigma_1^+I_k -S),$$ and therefore, \begin{equation}\label{equ6.5} \sigma_k(S) = \sigma_1^+ - \sigma_1(\sigma_1^+ I_k -S).\end{equation}Apply proposition~\ref{pro6.1} to $\sigma_1^+ I_k -S$ to obtain the  bounds on the approximation of the singular value $\sigma_1(\sigma_1^+ I_k -S)$ by $\vec{v}^T(\sigma_1^+ I_k -S)^q \vec{v}$. Combine these bounds with (\ref{equ6.4}) and (\ref{equ6.5}) and obtain that $\sigma_k^+ \geq \sigma_k(S)$ and furthermore, that $\sigma_k^+ <\sigma_k(S) -t$ with a probability at most $p_{k, \Theta, q}$, for $t= (\sigma_1^+ -\sigma_k)\Theta$. By letting $\Theta = t/ (\sigma_1^+ -\sigma_k)$, observing that $\sigma_k \geq 0$, so that $(\sigma_1^+  \geq \sigma_1^+ -\sigma_k$, applying  (\ref{equ6.3}), we verify that $\sigma_k^+$ satisfies both assumed requirements a) and b). \hfill $\Box$ The cost of performing of algorithm~\ref{alg6.1} is bounded by $(M(S) + O(k))q$, not counting the cost of the estimation of $\sigma_1^+$ and the generation of a random vector $\vec{v}$. Since the cost is proportional to $q$ and therefore, to $\sigma_1^+$, it could pay off if one may improve the bound $\sigma_1^+$ by performing the auxiliary computations based on proposition~\ref{pro6.1} for the same vector $\vec{v}$ and reusing the computed products $A^q\vec{v}$ at stage 4 of the algorithm. \section{Some properties of structuredmatrices}\label{sec7} \setcounter{equation}{0}  In the following, we will deal with various structured matrices, in particular, with Toeplitz and Hankel matrices, $$ T=(t_{i,j}), \  \  t_{i+1,j+1}=t_{i,j}, \  \   {\rm for \  all} \ i, j, $$$$ H=(h_{i,j}), \  \  h_{i+1,j-1}=h_{i,j}, \  \   {\rm for \  all} \ i, j, $$whose entries are invariant in their shifts into the diagonal and antidiagonal directions, respectively. These two classesof matrices are closely related to each other. In particular, $HJ$and $JH$ are Toeplitz matrices, and $JT$ and $JT$ are Hankel matrices,where $J$ is the reflection Hankel matrix, $J=\left( \begin{array}{cccc}0  & .  & 0  & 1 \\.  & .  & 1  & 0 \\.  & .  & .  &  \\0  & 1  & .  & . \\1  & 0  &    & 0 \end{array} \right)$, $J[v_0, \ldots, v_n]^T=[v_n, \ldots, v_0]^T$. An $N \times N$ Toeplitz or Hankel matrix is completely defined by the$2N-1$ entries of its two rows or two columns (the first and the last) andeven by a single row or column if the matrix is triangular. Clearly, a Hankel matrix is symmetric, and the transpose of a Toeplitzmatrix is a Toeplitz matrix. We immediately verify the next result. \begin{proposition}\label{pro7.1} For two polynomials $u(x)$ and $v(x)$ of (\ref{equ2.2}), their product $p(x)=u(x)v(x)$ has the coefficient vector $\vec{p}=U\vec{v}$, where $U=U(\vec{u})$ is the $(m+n+1) \times (n+1)$ Toeplitz matrixwith the first row $[u_0, 0, \ldots, 0]$ and the first column $[u_0, \ldots, u_m, 0, \ldots, 0]^T$.\end{proposition} \begin{corollary}\label{cor7.1} The multiplication of an $N \times N$ Toeplitz matrix $T$ by a vectorcan be reduced to multiplication of two polynomials of degrees $O(N)$and performed in $O(N\log N)$ ops. Furthermore,the latter cost bound also applies to multiplication by avector of $N \times N$ matrices of related classes of structured matrices, including Hankel matrices, $H=JT$, and the products of $O(1)$ Toeplitz or Hankel matrices or block matriceswith $O(1)$ blocks, which are all Toeplitz and/or Hankel matrices. \end{corollary} The latter class includes $(m+n) \times (m+n)$ {\it Sylvester(resultant) matrices} $S=S(\vec{u}, \vec{v})=[U(\vec{u}), V(\vec{v})]$, that is, $1\times 2$ block matrices formed with the pair of Toeplitz blocks,  $U(\vec{u})$ and $V(\vec{v})$ defined as in proposition~\ref{pro4.1}, so that thecoefficient vector of the polynomial $\beta(x) u(x) + \gamma(x) v(x)$ is given by$S\left[ \begin{array} {c} \vec{\beta}\\ \vec{\gamma} \end{array} \right]$. We will deal also with the $(m+n-k) \times (m+n-2k)$ submatrices$S^*_k$ and the $(m+n-2k) \times (m+n-2k)$ submatrices$S_k$ of $S=S_0$, called {\it subresultant matrices}, obtained bydeleting the last $k$ and $2k$ rows of matrix $S$, respectively, aswell as the last$k$ columns of each of its blocks  $U(\vec{u})$ and $V(\vec{v})$, for $k=0, \ldots, m$. Besides multiplication by a vector, other matrix operations are alsodramatically accelerated for structured matrices (cf. \cite{bp94}). In particular we have \begin{proposition}\label{pro4.2} For two polynomials $u(x)$ and $v(x)$ of (\ref{equ2.2}), their product $p(x)=u(x)v(x)$ has the coefficient vector $\vec{p}=U\vec{v}$, where $U=U(\vec{u})$ is the $(m+n+1) \times (n+1)$ Toeplitz matrixwith the first row $[u_0, 0, \ldots, 0]$ and the first column $[u_0, \ldots, u_m, 0, \ldots, 0]^T$.\end{proposition} 