From:	IN%"sad@cunyvms1.gc.cuny.edu" 23-MAR-1995 14:37:02.89To:	IN%"vpan@lcvax.lehman.cuny.edu"CC:	Subj:	This is the deterministic polynomial zeros.Return-path: <sad@cunyvms1.gc.cuny.edu>Received: from vice1.lehman.cuny.edu by lcvax.lehman.cuny.edu (PMDF #12809) id <01HOH8ZK74BK000O1W@lcvax.lehman.cuny.edu>; Thu, 23 Mar 1995 14:36 ESTReceived: from CUNYVMS1.GC.CUNY.EDU by vice1.lehman.cuny.edu (5.51/SMI-3.0DEV3) id AA05852; Thu, 23 Mar 95 14:27:07 ESTReceived: by cunyvms1.gc.cuny.edu (MX V3.3 VAX) id 31935; Thu, 23 Mar 1995 14:26:14 ESTDate: Thu, 23 Mar 1995 14:25:58 ESTFrom: sad@cunyvms1.gc.cuny.eduSubject: This is the deterministic polynomial zeros.To: vpan@lcvax.lehman.cuny.eduMessage-id: <0098DCB2.7B3748E6.31935@cunyvms1.gc.cuny.edu>\documentstyle[12pt]{article}\topmargin-.25in\oddsidemargin.05in\evensidemargin.05in\textheight8.55in\textwidth6in\newtheorem {algorithm}{\sc \bf Algorithm}[section]\newtheorem {problem}{\sc \bf Problem}[section]\newtheorem {theorem}{\sc \bf Theorem}[section]\newtheorem {proposition}{\sc \bf Proposition}[section]\newtheorem {lemma}{\sc \bf Lemma}[section]\newtheorem {corollary}{\sc \bf Corollary}[section]\newtheorem {fact}{\sc \bf Fact}[section]\newtheorem {remark}{\sc \bf Remark}[section]\newtheorem {definition}{{\sc \bf Definition}}[section]\newtheorem {claim}{\sc \bf Claim}[section]\def\theequation{\arabic{section}.\arabic{equation}}\begin{document}\begin{center}{\bf DETERMINISTIC IMPROVEMENT OF COMPLEX POLYNOMIAL FACTORIZATIONBASED ON THE PROPERTIES OF THE ASSOCIATED RESULTANT}\end{center} \begin{center}{\bf Victor Y. Pan}\\Mathematics  and Computer Science  Department \\Lehman College, CUNY \\Bronx, NY 10468 \\VPAN@LCVAX.BITNET \\(Supported by NSF Grant CCR 9020620 and PSC CUNY Award 664334)\end{center}\vspace{.5in} \noindent{\bf Summary.}We show an improvement of the known best deterministic approach tofactorization of a univariate polynomial $p(x)$ over the complex field.Previously, such an improvement was achieved by means of a randomperturbation of the coefficients of the input polynomial. We now applysome new techniques, which lead to a similar result without anyrandomization. The approach is unusual in exploiting some properties of theresultant of $p(x)$ and $p'(x)$, which is completely distinct from thetechniques usually appliedin order to computecomplex polynomial factorization.\vspace{.1in} \noindent{\bf Key Words:} complex polynomial factorization, resultant,derandomization\vspace{.3in} \noindent{\bf 1991 Math. Subject Classification:} 13P05, 68Q25, 68Q40, 65H05, 65Y20\newpage\noindent{\bf 1. Introduction} We consider the problem of approximate factorization of a univariatepolynomial of degree $n$, in the complex field,$$p(x) =\sum^n_{i=0} p_i x^i =  p_n \prod^n_{j=1} (x-z_j)\ , \ \ \ \ \parallel p(x)\parallel =\sum_i |p_i| \ , \eqno (1.1)$$that is, the problem ofcomputing approximations $z_j^*$ to the zeros $z_j$ of$p(x)$ such that$$\parallel  p(x) - p_n \prod^n_{j=1} (x-z^*_j) \parallel \leq\epsilon \ ,\eqno (1.2)$$for a fixed positive $\epsilon$. Here and hereafter, we assume that$\parallel \displaystyle{\sum_i}u_i x^i \parallel = \sum_i |u_i |$. Our initial idea was aboutperturbing slightly the coefficients of $p(x)$ so asto regularize the disposition of the zeros of $p(x)$ and to facilitatetheir subsequent approximation. Indeed, some of the most effective techniquesfor complex polynomial factorization via the isolation of the zeros ortheir clusters from each other \cite{p87}, \cite{p94a},\cite{r}, \cite{sc82} slow down their progress when thezeros are clustered in a certain irregular way, so that the isolationalgorithms peel them out one by one, rather than  isolate a considerablefraction of allthe zeros in one recursive step. A similar problem arises inthe QUICKSORT algorithm, but the analysis shows that the slowdown ofsorting occursrarely, and in the average case the algorithm runs fast. The analysis in\cite{p94} shows that a somewhatsimilar phenomenon takes place in the case offactorization of a complex polynomial, that is, the isolation algorithms of\cite{p87} and \cite{p94a} enable us to isolate (with a high probability) aconsiderable fraction of all the zeros already in the first recursive stepif we randomly perturb the $x$-free term of $p(x)$ (compare our Remark A.2in Appendix A). In the present paper we apply some new technique, which enables us toderandomize the algorithm of \cite{p94}, without affecting the asymptoticestimates for the parallel arithmetic complexity of \cite{p94}. Thus,we made the results of [P94] deterministic and improvedthe factorization algorithms of \cite{p87},\cite{p94a} and [Sc]by ensuring a faster isolation of the zeros of$p(x)$.Specifically, like [P94], we exploit the correlations, on the one hand,between the fast isolation of the zeros of $p(x)$ by means of application of Weyl's(quadtree) geometric construction and the magnitude of the $m$-thdiameter of the set $Z$ of the zeros of $p(x)$ (where, say, $m <n/2$),that is, the minimum diameter of all the subsets of cardinality $m$ ofthe set $Z$ [for formal definition see (2.2) in our section 2], andon the other hand, between the $m$-th diameter and the magnitude of theresultant $R=R(p(x),p' (x))$ of $p(x)$ and $p' (x)$. The main result of[P94] then followed from a probabilistic lower bound on the magnitude of$R(p(x)+\Delta ,p' (x))$, for a small random $\Delta$. (Thus $\Delta$defines a small random perturbation of $R$ via a small random change ofthe $x$-free term of $p(x)$.) Our present derandomization relies onobtaining deterministic lower bounds on $|R|$ after a similar butdeterministic perturbation (Corollary 3.2). As in [P94], the acceleration can be further accentuated for a largeclass of polynomials $p(x)$ if we also allow to perturb other coefficientsof $p(x)$. The influence of this perturbation on the $m$-th diameter isexpressed in terms of the partial derivatives of $R$ with respect to thecoefficients of $p(x)$, and in the present paper we established some newcorrelations between $R$ and its partial derivatives, which should motivatethe application of this heuristic extension of our main algorithm.[Compare Lemma 2.4, Remarks 2.1 and 6.1, Corollary 3.2 and the equations of(3.3) and recall thatcomputation of $R$ can be immediately extended to computation of allpartial derivatives of $R$ (see [L] or [BS]).] The proposed techniques work by accelerating the process of recursivelysplitting $p(x)$ into smaller degree factors as soon as the zeros of $p(x)$are partitioned into two or several groups isolated from each other.The basic splitting techniques can be traced back at least to [DH] andare currently best presented in the seminal paper [Sc82], which, however,is still available only as a very long unpublished manuscript with someproofs and many details missing. Thus, to make our paper more complete,we reproduce some relevant results from [Sc82] in Appendices A and B;we also supply some proofs and details missing in [Sc82] and simplify someblocks of the splitting algorithm. This part can be of some independentinterest for the study of the techniques for splitting a complexpolynomial into factors; the importance of such splitting has beenwell recognized. Summarizing, our study contributes to the four important areas: ofderandomization, study of the resultants, of their correlation to fastpolynomial factorization in the complex field and the techniques forsplitting a complex polynomial into factors, where the zeros ofeach factor are isolated from the other zeros of $p(x)$. Our main result is stated in the next theorem.[Hereafter, $O_A (t,p)$ shows that the computation can be performedusing time $O(t)$ and simultaneously $O(p)$ processorson EREW arithmetic PRAM;we assume thatthis also implies $O_A (tp,1)$, that is, time bound $O(tp)$ with a singleprocessor (Brent's principle).] \noindent{\bf Theorem 1.1}{\it Given the complex or real coefficients$p_0 ,p_1 , \ldots , p_n$, $p_n \neq 0$, of a polynomial $p(x)$ of (1.1),and a positive $\epsilon$satisfying the relations}$$\parallel p(x) \parallel = 1 \ , \ \ \ 1/4 \leq r(C,p(x)) = \max_j|z_j -C| \leq 1 \ ,$$$$C = (1/n) \sum_j z_j =-p_{n-1} /(np_n ) \ ,$$$$\epsilon < n|p_n |20^{-n} \ , \eqno (1.3)$${\it it is possible to compute the factorization of $p(x)$ (withinthe error norm at most $\epsilon$) into the factors of degrees less than$n/2$ at the asymptotic parallel cost$$O_A ( \phi (n, \epsilon ) \log n ,  \ n^3 / \phi (n, \epsilon )) \ ,$$where$$\phi (n, \epsilon ) =\log ( \log (1/ \epsilon ) +(n/ \epsilon^{1/n}) ) \ .$$}\vspace{.1in} \noindent{\bf Remark 1.1}The inequality (1.3) has been based on a result from [Sc82] that related (1.2) to a certain upper bound on the errors $|z_j^*-z_j|$ of approximation to the zeros $z_j$ of $p(x)$by $z_j^*$ for $j=1,2, \ldots , n$. As soon as one deduced from (1.2) a smaller bound on these errors, we could soften or relax the bound (1.3) on$\epsilon$. \noindent {\bf Remark 1.2.} We refer the reader to the recent papers  [N], [NR] and [P], on an alternative balanced splitting approach to factorization of complex univariate polynomials.\vspace{.2in} \noindent{\bf2. Some Auxiliary Properties of the Resultant of a Polynomial andIts Derivative} The resultant $R=R(p,p')$ of a polynomial $p(x)$ of (1.1) and of itsderivative $p'(x)$ is defined \cite{vdW} as the determinant of the associated$(2n-1) \times (2n-1)$ Sylvester matrix,$$R ={\rm det} \  S  , \ S =\left[\begin{array}{cccccccc}p_n &~ &O &~ &(np_n) &~&~ &O \\\vdots &\ddots &~ &\vdots & &\ddots &~ \\p_2 &\cdots &p_n &~ &(2p_2) &~ &~ &~\\p_1 &\cdots  &p_{n-1} &~ &p_1 &~ &\cdots &(np_n) \\p_0 &\cdots &p_{n-2} &~ &~ &~ &~&~ \\~ &\ddots &\vdots &~ &~ &~&\ddots &\vdots \\O &~ &p_0 &~ &O &~ &~&p_1\end{array}\right] \ .\eqno (2.1)$$ We will next observe some simple properties of the resultant $R$. {\bf Observation 2.1} $R$ {\it is a polynomial in $p_0 , p_1 , \ldots ,p_n$ with integer coefficients. It has degree $n-1$ in $p_0$, hasdegrees at most $n$ in each of $p_1 , \ldots , p_n$, and has totaldegree $2n-1$.} We will only cite the next two lemmas in Remark 6.1, but theselemmas themselves may be of independent interest. \noindent{\bf Lemma 2.2} \cite{vdW}. $R=(-1)^{(n-1)n/2} p_n^{2n-1} \displaystyle{\prod_{1 \leq i <j \leq n} (z_i -z_j )^2}$. \noindent{\bf Lemma 2.3}{\it Let $R_t$ denote $R(g,g')$, where $g(x) =p(tx) = \displaystyle{ \sum_{i=0}^n} g_i x^i$, $g_i = p_i t^i$. Then$R_t = t^{a(n) }R$, $a(n)=(2n-1)n-(n-1)n =n^2$.} {\em Proof.}Lemma 2.3 follows from Lemma 2.2 since $g(x)$ has zeros$z_i /t$, $i=1, \ldots ,n$. $\Box $ \noindent{\bf Lemma 2.4}.$\displaystyle{ R = {1 \over n^2} \sum_{k=1}^n (kp_k ){\partial R \over \partial p_k} }$. {\em Proof.}Recall that $\displaystyle{ {dR_t \over d t}= \sum_{k=0}^n \left( {\partial R_t \over \partial g_k} \right)\left( {dg_k \over dt} \right)}$, set $t=1$, and verify that in this casewe have $g_k =p_k$, ${\partial R_t \over \partial g_k} = {\partial R \over\partial p_k}$, and ${dg_k \over dt} =kp_k$ for all $k$. Combine allthese equations with Lemma~2.3 and arrive at Lemma 2.4. $\Box$ \noindent{\bf Remark 2.1}For $n=2$, $R=p_2 (4p_0 p_2 -p_1^2 )$, so that$R=0$ if $p_1^2 =4p_0 p_2$, whereas in the case, where $p_1^2 =4p_0 p_2$, all the first order partial derivatives,$\partial R/ \partial p_0 =4p_2^2$, $\partial R /\partial p_1 =-2p_1 p_2$ and $\partial R/\partial p_2 =8p_0 p_2-p_1^2$, are nonzeros unless $p_0 p_1 p_2 = 0$. In particular, if we shift from $p(x)$ to $q(y) = p(y+a) = q_2 y^2 +q_1 y+q_0$, for a sufficiently large $a$, we would not change $R$ but would ensure that all the partial derivatives of$R$ in $q_i$ do not vanish. Similar observations apply in the case where $n=3$.\vspace{.1in} Hereafter $d_m$ denotes the $m$-th diameter of the set $\{ z_1 , \ldots ,z_n \}$ defined as follows:$$d_m = \min_{|T|=m} \ \max_{i,j \in T} |z_i -z_j | \ , \ \ \ \m=2,3, \ldots , n \ , \eqno (2.2)$$where the minimum is over all the subsets $T$ of the set $\{ 1,2 , \ldots ,n \}$ having cardinality $m$, so that, in particular,$\displaystyle{d_n  = \max_{1 \leq i,j \leq n} |z_i -z_j |}$. By applying Lemma 2.2, we immediately relate the value$$D =D (p_0 , \ldots , p_n ) = |R/p_n^{2n-1} | \eqno (2.3)$$to the diameters $d_m$ and $d_n$ for $2 \leq m \leq n$. \noindent{\bf Lemma 2.5}{\it Let $ R \neq 0$, $2 \leq m \leq n$. Then $D \leqd_n^{(n-1)n} (d_m /d_n )^{(m-1)m}$}. By applying Lemma 2.5 for $n=2m$ and for $m= \lfloor n^a \rfloor$,$0 < a<1$,we obtain that$$d_{n/2} \geq D^{4/((n-2)n)} / d_n^{(3n-2)/(n-2)} \ ,\  \ {\rm for \ even} \ n\ , \eqno (2.4)$$$$d_{\lfloor n^a \rfloor} \geq\left( D /d_n^{(n-1)n-\lfloor n^a -1 \rfloor \lfloor n^a \rfloor}\right)^{1/(\lfloor n^a -1 \rfloor \lfloor n^a \rfloor )}\ , \ \  \ \ {\rm for} \ 0 < a < 1 \ . \eqno (2.5)$$\vspace{.3in} \noindent{\bf3.  Lower Bounds on the Magnitude of the Resultant} Later on we will fix an integer $h$, a positive $\rho$ and a complex$\Delta$ such that $0 \leq h \leq n$, $| \Delta | = \rho$. Hereafter wewill set$$p_h (\Delta ,x)=p(x)+x^h \Delta =p_{h,n} \prod_{k=1}^n (x-z_{h,k} (\Delta )) = \sum_{i=0}^n p_{h,i} (\Delta )x^i \ , \eqno (3.1)$$where $p_{h,h} =p_h + \Delta $, $p_{h,i} =p_i$ for $i \neq h$. Wewill extend Lemmas 2.2--2.5 and the definitions (2.1)--(2.3) byreplacing $p(x)$, $p_i$, $z_k$, $R$, $R_t$, $d_m$ and $D$ by$p_h (\Delta ,x)$, $p_{h,i} (\Delta ) ,z_{h,k} (\Delta )$,$R_h ( \Delta )$, $R_{t,h} ( \Delta )$, $d_{h,m} ( \Delta )$ and$D_h ( \Delta )$, respectively, and we will also represent $R_h ( \Delta )$as a polynomial in $\Delta$,$$R_h(\Delta ) = \sum^n_{i=0} r_{h,i}  \Delta^i\ , \eqno (3.2)$$where $r_{h,i}$ are polynomials in $p_0,\ldots,p_n$ with integercoefficients that have total degrees at most $2n-1-i$, $i=0,1,\ldots,n$.Furthermore, since $R_h(\Delta)$ is a polynomial in $p_h+\Delta, \  p_0, \  p_1, \ldots , \  p_{h-1}, \ p_{h+1}, \ \ldots, p_n$,  we obtain that $$r_{h,1} = \partial R/\partial p_h \ , \qquad h = 0,1,\ldots,n \ .\eqno (3.3)$$Moreover, it is easily deduced from (2.1) that $$r_{0,n}=0 \ , \quad r_{0,n-1}  = (-1)^{n(n-1)/2} (n p_n)^n \ .\eqno (3.4)$$ Next, we will estimate $|R_h (\Delta )|$ for $\Delta$ ranging on the circle$\{ \Delta : | \Delta | = \rho \}$. \noindent{\bf Lemma 3.1} $\displaystyle{  \int_0^{2\pi} |R_h( \rho \exp (2 \pi \ \sqrt{-1} \ \phi)) |d \phi \geq |r_{h,i}|\rho^i}$ {\it for}$i=0,1, \ldots ,n$. {\bf Proof} \cite{p94}. It is well-known from the analytic function theory that$$\int_{|\Delta |=\rho} (R_h (\Delta )/\Delta^{j+1} )d\Delta = 2\pi\ \sqrt{-1} \ r_{h,j} \ .$$By representing $\Delta $ as $\rho\exp (2\pi \ \sqrt{-1} \ \phi )$, weobtain that$$(1/ \rho^j ) \int_0^{2 \pi} |R_h (\Delta )|d \phi = \int_0^{2 \pi}|R_h (\Delta )/\Delta^{j+1} |\rho \ d\phi \geq 2\pi |r_{h,j}|,$$which proves Lemma 3.1. $\Box$ We now observe that $|R_h (\Delta )| \leq \sum^n_{i=0} |r_{h,i} \rho^i |\leq (n+1) \max \{ |r_{h,i} |\rho^i \}$, combine this upper boundon $|R_h ( \Delta )|$ with the bound of Lemma 3.1, and obtain the followingresult. \noindent{\bf Corollary 3.2}{\it If for a fixed positive $\rho$, a random $\Delta$is chosen on the circle $\{ \Delta :|\Delta | = \rho\}$ under theuniform probability distribution, then$$|R_h (\Delta )| > 0.5 |r_{h,i}| \rho^i$$with a probability at least $1/(2n+1)$, for any pair of integers $i$ and$h$ fixed in the range from $0$ to $n$.} Randomized versions of ourfactorization algorithms may rely on the bound of Corollary 3.2,as in [P94].We will  present their derandomized versions based on the followingresult. \noindent{\bf Lemma 3.3}{\it Let $M$ and $\theta$ denote two fixed positiveconstants, let $\displaystyle{ P(y) =\sum_{i=0}^n P_i y^i}$ denote apolynomial,  let $|P(y)|>M$ on an unknown set $Y$of points on the unitcircle $\{ y: |y|=1 \}$, and let this set $Y$ have Lebesque's measure atleast $\theta$. Then$$\max_{0 \leq g < N} |P(y_g )|>M/ \sqrt{2} \ ,$$where $y_g = \exp (2 \pi \ g \ \sqrt{-1} / N)$; $N> \pi n/ \theta$,$g=0,1, \ldots ,N-1$}. {\em Proof.}With no loss of generality, we may assume that $|{\rm Re} \P(y)|>M/ \sqrt{2}$ on a subset $Y_{Re}$ of $Y$ having Lebesgue's measure atleast $\theta /2$. [Otherwise$|{\rm Im} \ P(y)|=|{\rm Re} (\sqrt{-1} \ P(y))|$ has this property,and we may shift from $P(y)$ to $\sqrt{-1}$ $P(y)$.] Since $P(y)$ isa polynomial of degree at most $n$, the set $Y_{Re}$ is the union of at most$2n $ arcs on the circle $\{ y: \ |y|=1 \}$.(To verify this property, substitute$y = t+ \sqrt{t^2 -1}$ and $y=t- \sqrt{t^2 -1}$ and represent $P(y)$ inthe form of $P_0 (t) + \sqrt{-1} \ P_1 (t)$ for $-1 \leq t \leq 1$,$|y| =1$, where $P_0 (t)$ and $P_1 (t)$ are polynomials of degrees atmost $n$.) The longest of these arcs has length at least $\theta /(4n)$.Therefore, it includes at least one point $y_g$, since the arc between $y_g$and $y_{g+1}$ has length $2\pi /N \leq \theta /(4n)$, for every $g$.$\Box$ Combining Corollary 3.2 and Lemma 3.3, for $\Delta =\rho y$, $P(y)=R( \rho y)$,$M=0.5 |r_{h,i}|\rho^i$, $\theta =1/(2n+1)$, gives the next result. \noindent{\bf Corollary 3.4}$\displaystyle{\max_{0 \leq g <N}|R_h (\rho y_g )|> 0.5 |r_{h,i} |\rho^i / \sqrt{2} }$ {\it for any pair ofintegers $h$ and $i$,lying in the range from 0 to$n$, for $y_g = \exp (2\pi g \ \sqrt{-1} /N)$, andfor any integer $N \geq \pi n(2n+1)$.} \noindent{\bf4. Polynomial Factorization via Isolation of Zeros} We will first recall some definitions from \cite{p87}, \cite{p94a}. Hereafter,$D(c,r)$ and $S(c,r)$ denote the closed disc and square, respectively,lying on the complex plane and having the center $c$ and either radius$r$ or the side length $2r$, respectively. For $ R \geq r$ the value$R/r$ is said to denote the {\it isolation ratio of} $D (c,r)$ or$S(c,r)$ if $R$ equals the maximum positive number such that the disc$D(c,R)$ (respectively, the square $S(c,R)$) contains exactly as manyzeros as $D(c,r)$ (respectively, $S(c,r)$), (compare Figure 1). Discs andsquares are called $f$-{\it isolated} if their isolation ratios are atleast $f$. {\it Weyl's algorithm} [He](also known as the {\it quadtree algorithm})approximates allthe zeros of a polynomial $p(x)$ lying in  a fixed suspect square $S$on the complex plane. (One may start with asquare $S$ containing eitherall the zeros of $p(x)$ or a part of them.) Thealgorithm consists of a sequence of recursive steps. Every step startswith a set of suspect squares, partitions each of them into 4congruent subsquares, and approximates the distance from their centersto the nearest zeros of $p(x)$. If this proximity test proves that asubsquare contains no zeros of $p(x)$, then this subsquare isdiscarded; otherwise it is called suspect and is further partitionedat the next recursive step (see Figures 2 and 3). The input square $S=S(C,r)$ containing all the zeros of $p(x)$ can bedefined by its readily available center,$$C ={ 1 \over n} \sum_{j=1}^n z_j =p_{n-1} /(np_n ) \ , \eqno (4.1)$$and by the value $r$, computed at the cost $O_A (\log n,n)$ of performing $O(1) FFT$s, at $O(n)$ points each,\cite{tu}, and satisfying the relations$$\theta r \leq r(C,p(x)) = \max_j |C-z_j | \leq r \ , \eqno (4.2)$$where $\theta$ is a fixed constant, $0 < \theta < 1$, and $r(C,p(x))$ denotes the distance from $C$ to the farthest zero of $p(x)$. By taking thisinto account and by using Weyl's construction, together with some customarytechniques, such as winding number computation via FFT and Graeffe'smethod for root squaring, we obtain the next result. \noindent{\bf Lemma 4.1}{\it Let a complex $C$ and positive $r$ and $\theta$satisfy the above relations (4.1) and (4.2). Then for some $H=H(n,r,\theta ) = O(\log n)$ and for every natural $K>H$,  one may, at the cost$O_A (K \log n ,n^2 )$, compute complex $C_i$, positive $r_i$, andnatural $g \geq 2$ and $k_i \geq 1$, $i=1, \ldots ,g$, such that for all$i$ we have$$r_i \leq rn/2^K$$and$S(C_i ,r_i )$ is 8-isolated and contains exactly $k_i$ zeros of $p(x)$;furthermore, $\displaystyle{ \sum_{i=1}^g k_i =n}$, and the set$S( C_j ,r_j ) \cap S(C_h ,r_h )$ is empty for all pairs of $j$ and $h$such that $1 \leq j <h \leq g$.} Having the squares $S(C_i ,r_i)$ of Lemma 4.1 available, we may easilycompute the associated factorization of the polynomial $p(x)$. \noindent{\bf Lemma 4.2}{\it Let $S(C_i ,r_i )$ denote the squares defined byLemma 4.1 for $i=1, \ldots ,g;g \leq n$. Let these squares be availableas well as the associated integers $k_i$ and a positive $\epsilon^*$.Let $b =b( \epsilon^* ) = \log (1/ \epsilon^* )$. Then at the cost$O_A ( \log n \log (bn)$, $gn \log b/\log (bn))$, one may compute thecoefficients of $g$ monic polynomials $p_1 (x) , \ldots ,p_g (x)$ suchthat $p(x)$ has exactly $k_i$ zeros, all lying in $S(C_i ,r_i )$,$i=1, \ldots ,g$, and }$$\parallel p (x) -p_1 (x) \ldots p_g (x) \parallel \leq \epsilon^* \ .\eqno (4.3)$$  Lemma 4.2 is deduced by applying thealgorithms of \cite{sc82}modified in Appendices A and B to the discs $D(C_i ,r_i \sqrt{2})$; thedisc $D(C_i ,r_i \sqrt{2})$ contains exactly $k_i$ zeros of $p(x)$ andis $(3/ \sqrt{2})$-isolated. To obtain factorization of $p(x)$ into linear factors, one may apply thefactorization algorithm supporting Lemma 4.2, at first to $p(x)$ and thenrecursively to $p_i (x)$ for $i=1, \ldots ,g$. The inspection reveals thatthis recursive process is the slowest in terms of both numbers ofrecursive steps and arithmetic operations involved if $g=2$ and $k_1 =1$in all the recursive steps. On the other hand, by transition from $p(x)$ to a neighboring polynomial$p(x)+ x^i \Delta $for a small $|\Delta |$, we may ensure balanced partitionof the zeros of $p(x)$, so that $k_i < m$ for a fixed $m \leq n/2$ or$m \leq \lfloor n^a \rfloor$, $a<1$, say. We may apply the same technique tothe resulting approximate factors of $p(x)$, although the approximationgenerally deteriorates as the degrees of the factors decrease.\vspace{.3in} \noindent{\bf 5. Balancing the Splitting of the Set of the Zeros} We will next obtain the relations for $\rho$ and $K$ for which Weyl'salgorithm partitions the zeros of $p(x)$ into isolated groups with less than$n/2$ zeros in each group, thus balancing the degrees of the factors inthe subsequent splitting of $p(x)$. By combining the estimates of Lemmas 2.4 and 4.1, we obtain the followingresult. \noindent{\bf Lemma 5.1}{\it Let the integer $K$ and the side length $2r_i$ of thesquares $S(C_i ,r_i )$ of Lemma 4.1 satisfy the relation$$r_i \leq 2 \sqrt{2} rn /2^K < (D/d_n^{(n-1)n-(m-1)m})^{{1 \over(m-1)m}}\leq d_m$$for some natural $m <n$ and for $r$ of (4.2). Then every such a squarecontains less than $m$ zeros of $p(x)$.} We will apply Lemma 5.1 to the polynomial $p(\Delta ,x)$ of (3.1)replacing $D$ by $D_0 (\Delta )$, $d_j$ by $d_j ( \Delta )$ for $j=m$and $j = n$,and $r$ by $r_0 (\Delta )$, so that the assumption of Lemma 5.1 is rewrittenas follows:$$r_i \leq 2 \sqrt{2} \ r_0 (\Delta ) n /2^K < (D_0 (\Delta ) /d_n (\Delta)^{(n-1)n-(m-1)m)})^{{1 \over (m-1)m}} \leq d_m (\Delta ) \ .$$We need to choose $K$ and $\rho =|\Delta |$ so as to satisfy themiddle inequality, which we rewrite as follows:$$D_0 (\Delta )>(2 \sqrt{2} \ r_0( \Delta )n2^{-K})^{(m-1)m}  d_n (\Delta)^{(n-1)n-(m-1)m} \ , \eqno (5.1)$$where $d_n (\Delta ) \leq 2r_0 (\Delta )$. Now, by choosing $m=n/2$(assuming even $n$), so that $(m-1)m=(n-2)n/4$, and applying Corollary 3.2 for $h=0$, $i=n-1$ and the bound (3.4), we deducethat (5.1) holds with a high probability for a random choice of $\Delta$on the circle $|\Delta |=\rho$, provided that$$0.5|r_{0,n-1} /p_n^{2n-1}|\rho^{n-1} \geq d_n^{(n-1)n}(\Delta) n^{n(n-2)/4}2^L \eqno (5.2)$$for $L=(0.5-K)(n-2)n/4$ and for $r_{0,n-1}$ of (3.4). Alternatively, we may rely on Corollary 3.4, rather than 3.2, andensure (5.1) deterministically for$$L=(0.5-K)(n-2)n/4+0.5 \ , \ \ \ \\Delta =\rho\exp (2\pi g \ \sqrt{-1} /N ) \ , \eqno (5.3)$$ $N \geq \pi n(2n+1)$ and forat least one integer $g$, $0 \leq g <N$.\vspace{.3in} \noindent{\bf 6. Factorization Algorithm and a Discussion} By relying on (5.2), (5.3) and Lemmas 4.1, 4.2 and 5.1, we arrive at thefollowing algorithm. {\bf Algorithm 6.1.} {\bf Input:} coefficients $p_0 , \ldots ,p_n$of $p(x)$ and a positive $\rho$ satisfying the bounds$$r(C,p(x)) \leq 2r(C,p(x)+\Delta ) \leq 3 r (C,p(x)) \ {\rm if } \|\Delta | = \rho \ , \ \ C = -p_{n-1}/(np_n ) \ . \eqno (6.1)$$ {\bf Output:} a complex $\Delta$, $|\Delta |=\rho$, an integer $g \geq 2$and the coefficients of $g$ polynomials $p_i^* (x,\Delta )$ having degreesless than $n/2$ and satisfying the bound$$\parallel p(x) - \prod_{i=1}^g p_i^* (x, \Delta ) \parallel \leq 2 \rho \ .\eqno (6.2)$$ {\bf Computation.} \begin{enumerate}\item Evaluate $R_0 (\Delta )$, the resultant of $p(x)+\Delta$, for$\Delta = \rho\exp (2\pi j \ \sqrt{-1} /n)$, $j=0,1, \ldots ,n-1$.\itemApply scaling and inverse FFT(at $n$ points) to the values of $R_0 (\Delta )$, computed at stage 1, andevaluate the coefficients of $R_0 (\Delta )$ represented as a polynomial in$\Delta $ of degree at most $n-1$.\itemEvaluate $R_0(\Delta )$ for $\Delta = \rho\exp (2 \pi j \ \sqrt{-1}/N)$,$j=0,1, \ldots ,N-1$, for $N = \lceil \pi n (2n+1) \rceil$, by applyingforward FFT at $N$ points. Of these $N$ points $\Delta$, select one,$\Delta_0$, that maximizes $|R_0 (\Delta )|$.\item Apply $K$ recursive steps of Weyl's algorithm (supporting Lemma4.1) to the polynomial $p(x)+\Delta_0 $,where $K$ satisfies (5.2), (5.3).\item Apply the factorization algorithms of Appendices A and B in order tosplit $p(x)+\Delta_0$ into factors according to Lemma 4.2, applied (for$\epsilon^* =\rho )$ to the set of squares $S(C_i ,r_i )$ and integers $k_i$computed at the previous stage 4. Output the computed factors of$p(x)+\Delta_0$.\end{enumerate} The correctness of the algorithm follows from the lemmas and estimatesof the previous sections (see the details in \cite{p94}). The parallel complexityestimates supporting Theorem 1.1 also follow. [Apply the complexity estimatesof Lemmas 4.1 and 4.2, note that the evaluation of the resultant for each$\Delta$ at stage 1 can be done at the cost bounded by $O_A (\log^2 n,n^2 /\log n)$ \cite{p92} and recall that FFT at $s$ points costs$O_A (\log s,s)$.] \noindent{\bf Remark 6.1}The sequential time bound implied in the proof ofTheorem 1.1 can beimproved under the randomized model of computing(see [P94]). Furthermore,Lemma 5.1 shows that application to $p(x)$ of the algorithmssupporting Lemmas 4.1 and 4.2 should yield balanced recursive splitting of$p(x)$ for a large class of polynomials $p(x)$ (for which $D$ is not smalland $d_n$ is not large).[Note that $d_n$ is readily estimated and $D=|R/p_n^{2n-1}|$ is readilyavailable at stage 2 of Algorithm 6.1, since $R=R_0$ (0).]For many of such polynomials even splitting them into products of severalsmall degree polynomials can be achieved in a single splitting step[see (2.5)].Moreover, Lemma 2.3 and Remark 2.1 show that generally the value$\displaystyle{ \max_{1 \leq i \leq n} |\partial R /\partial p_i | }$is bounded away from 0 as soon as so does the value $R$. Therefore, we mayrecall (3.3) and Corollary 3.4 and deduce thatby shiftingto the polynomial $p(x)+x^i \Delta$ for an appropriate $i$ one may[for a large class of polynomials $p(x)$]improvethe error bound in Theorem 1.1, without affecting the asymptotic cost bound(see [P94]). \noindent{\bf Remark 6.2}The approach of this paper can be applied recursively(to every computed factor $p_i (x,\Delta ))$ in order to approximatethe lowest degree factors of $p(x)$. However, as the degrees of theinput polynomials decrease in each splitting step, the bounds on the errorsof approximation grow a little unless we allow to increase the number ofsteps $k$ in Weyl's algorithm [compare (5.2) assuming that $n$ decreases].  \normalsize\begin{thebibliography}{92} \footnotesize \bibitem[Ah]{ah}L. Ahlfors, {\bf Complex Analysis}, McGraw-Hill, New York, 1979.\bibitem[AHU]{ahu}A.V. Aho, J.E. Hopcroft, J.D. Ullman, {\bf The Design and Analysis ofComputer Algorithms}, Addison-Wesley, Reading, Mass., 1976.\bibitem[BP86]{bp86}D. Bini, V.Y. Pan, Polynomial Division and Its ComputationalComplexity, {\bf J. of Complexity}, 2, 179-203 (1986).\bibitem[BP94]{bp94}D. Bini, V.Y. Pan, {\bf Polynomial and Matrix Computations}, vol. 1, Birkhauser, Boston, 1994.\bibitem[BS]{bs}W. Baur, V. Strassen, On the Complexity of PartialDerivatives, {\bf Theor. Computer Science}, 22, 317-330, 1983.\bibitem[DH]{dh}B. Dejon, P. Henrici (editors), {\bf Constructive Aspects of theFundamental Theorem of Algebra}, Wiley, New York, 1967.\bibitem[He]{He}P. Henrici, {\bf Applied and Computational Complex Analysis}, Wiley,New York, 1974.\bibitem[L]{l}S. Linnainmaa, Taylor Expansion of the Accumulated Rounding Errors,{\bf BIT}, 16, 146-160, 1976.\bibitem[KP]{kp}E. Kaltofen, V.Y. Pan, Processor Efficient Parallel Solution of LinearSystems over an Abstract Field, {\bf Proc. 3rd Annual ACM Symp. onParallel Algorithms and Architectures}, 180-191, ACM Press, New York, 1991.\bibitem[M]{m}M. Mignotte, An Inequality about Factors of Polynomials, {\bf Math.of Computation}, 28, 1153-1157, 1974.\bibitem[N]{n}C.A. Neff, Specified Precision Polynomial Root Isolation Is in NC,{\bf Journal of Computer and System Sciences}, 48, 429-463, 1994.\bibitem[NR]{nr}C.A. Neff, J.H. Reif, An $O(n^{1+\epsilon} \log b)$ Algorithm for theComplex Root Problem, {\bf Proc. 35th Ann. IEEE Symp. on Foundations ofComputer Science}, 540--547, IEEE Computer Society Press, 1994.\bibitem[P87]{p87}V.Y. Pan,Sequential and Parallel Complexity of Approximate Evaluation of PolynomialZeros,{\bf Computers and Math. (with Applics.)}  14 (8)591--622 (1987).\bibitem[P90]{p90}V.Y. Pan, Parallel Least-Squares Solution of General and Toeplitz-likeLinear Systems, {\bf Proc. 2nd Ann. ACM Symp. on Parallel Algorithmsand Architectures}, 244-253, ACM Press, New York, 1990. \bibitem[P92]{p92}V.Y. Pan,Parametrization of Newton's Iteration for Computations with StructuredMatrices and Applications, {\bf Computers and Math. (with Applics.)},24 (3)  61--75 (1992).\bibitem[P94]{p94}V.Y. Pan,New Resultant Inequalities and Complex PolynomialFactorization, {\bf SIAM J. on Computing}, 23 (5) 934--950 (1994).\bibitem[P94a]{p94a} V.Y. Pan, New Techniques for Approximating ComplexPolynomial Zeros, {\bf Proc. 5th Ann. ACM-SIAM Symp. on DiscreteAlgorithms}, 260-270, ACM Press, New York, andSIAM Publications, Philadelphia, 1994.\bibitem[P]{p}V.Y. Pan,Optimal and Nearly Optimal Algorithms for Approximating Complex Polynomial Zeros, submitted to {\bf Computer and Math. (with Applics.).}\bibitem[R]{r}J. Renegar, On the Worst-case Arithmetic Complexity of ApproximatingZeros of Polynomials, {\bf J. of Complexity}, 3 (2) pp. 90-113 (1987).\bibitem[Sc82]{sc82}A. Sch\"{o}nhage, The Fundamental Theorem of Algebra in Terms ofComputational Complexity, Manuscript, {\bf Dept. of Math., Univ. ofT\"{u}bingen}, T\"{u}bingen, Germany, 1982.\bibitem[Sc85]{sc85}A. Sch\"{o}nhage, Quasi-GCD Computation, {\bf J. of Complexity}, 1,118-137 (1985).\bibitem[Tu]{tu}P. Turan, {\bf On a New Method of Analysis and Its Applications},Wiley and Sons, New Jersey, 1984.\bibitem[vdW]{vdW}B.L. van der Waerden, {\bf Modern Algebra}, Frederick UngarPublisher, New York, 1953.\end{thebibliography}\appendix\normalsize{\normalsize\bf Appendices}\def\theequation{\Alph{section}.\arabic{equation}}\normalsize \section{\normalsize\bfSplitting  a Polynomial into Factors}\label{appA}\setcounter{equation}{0}  In this appendix we will assume an initial isolation ratio of$f\geq 4$ (or we could  have similarly choose $f\geq 8$, say), forthe input squares of the algorithm supporting Lemma 4.2and will show how to split the polynomial $p(x)$into  factors,  according to (4.3) for$\varepsilon^*=2^{-b}$.Moreover, we will consider the extension of splitting to each of thefactors, assuming again that for each of them the assumptions of Lemma4.2 are satisfied, which can be achieved by means of Weyl's algorithmand other techniques cited in section 4. In this and the next appendiceswe will assume that the recursive splitting continues until we arrive atthe $n$ linear factors of $p(x)$ and will show that all splitting stages(supporting Lemma 4.2) can be performedat the overall sequential computational cost \begin{equation}\label{equaA.1}c_A= O_A (n^2\log n\log (bn)).\end{equation} \noindentBy using the algorithms of the next appendix, one may achieve processorefficient acceleration to reach parallel time $O((k+\log (bn)) \log n)$,where $k$ is the number of splitting steps required. In particular, for $k=1$, we arrive at the following result, which is not used in this paper but is of independent interest.\begin{proposition}\label{propA.1} Let for polynomial $p(x)$ of (1.1) a fixed disc $D$ be $f$-isolated for $f=4$ (or for any fixed constant $f>1$). Then the coefficients of the factors $F(x)$ and $G(x) = p(x)/F(x)$ of $p(x)$ having all their zeros inside and outside $D$, respectively, can be approximated within $2^{-b}$ at the cost $O_A(\log n \log (bn), n)$.\end{proposition} The splitting algorithm can be traced back at least to[DH], pp.~193--208. We will follow its more recentversion of \cite{sc82}, although our presentation andanalysis are a little simplified against \cite{sc82}; inparticular, it is sufficient for our purpose toperform splitting over an$f$-isolated disc for $f\geq 4$,versus $f\geq 1+1/n$ in  \cite{sc82}, and therewe do not need to use any advanced techniques  appliedin \cite{sc82} to the error analysis of approximatingthe power sums of the zeros of $p(x)$. With no loss of generality, we will assume that the$f$-isolated input disc (for $f\geq 4$) is $D=D(0,1/2)$and that this  disc as well as the disc $D(0,2)$contains exactly$k$ zeros of $p(x),\ z_1,\ldots, z_k,\ 1\leq k<n$.Furthermore, denote that\begin{equation}\label{equaA.2}p(x)=F(x)G(x),\ \  F(x)=\prod_{j=1}^k(x-z_j).\end{equation}Given $p(x)$, the evaluation of $F(x)$ and $G(x)$ is called{\em splitting $p(x)$ over the disc $D$}. By recursivelysplitting at first $p(x)$, then $F(x)$ and $G(x)$, and so on,we may finally arrive at linear factors of $p(x)$ or, moregenerally, at the factors having all their zeros in some$f$-isolated $\varepsilon_h$-neighborhoods (with$\varepsilon_h<\varepsilon$ for all $h$),which are then output as the solution to theproblem of approximating $z_1,\ldots, z_n$. According to a standard pattern for splitting,we approximate $F(x)$ and $G(x)$ in two stages.At first we apply the FFT based numericalintegration to compute relatively roughapproximations $F_0(x)$ to $F(x)$ and  $G_0(x)$to $G(x)$; then we apply Newton's iteration tocompute a sequence of recursively improvedapproximations $F_h(x)$ and $G_h(x)$,$h=1,2,\ldots $. The first stage of this splitting processbegins with the evaluation of\begin{equation}\label{equaA.3}s_g^*=\frac{1}{Q}\sum_{q=0}^{Q-1}\omega^{(g+1)q}\frac {p'(\omega^q)}{p(\omega^q)},\ g=1,2,\ldots, Q-1,\end{equation}where $\omega=\exp (2\pi\sqrt {-1}/Q)$ is aprimitive $Q$-th root of 1 andwhere $Q$ will be chosen later on. The values$s^*_g$ approximate the power sums$s_g^{(k)}=\sum_{i=1}^kz_i^g$ of the zeros of$p(x)$, lying inside the disc $D(0,1/2)$. Toestimate the approximation errors, consider theLaurent expansion:\begin{equation}\label{equaA.4}p'(x)/p(x)=\sum_{j=1}^n\frac{1}{x-z_j}=-\sum_{g=1}^{\infty}s_{-g}^{(k)}x^{g-1}+\sum_{g=0}^{\infty}s_g^{(k)}x^{-g-1}=\sum_{h=-\infty}^{\infty}c_hx^h.\end{equation}Here, $|x|=1$, $s_g^{(k)}=\sum_{j=1}^k z_j^g$;$s_{-g}^{(k)}=\sum_{j=k+1}^n(1/z_j^g)$.Consequently, $s_g^{(k)}$ is the $g$-th powersum of all the $k$ zeros of $p(x)$ lyinginside the disc $D(0,1/2)$, whereas $s_{-g}^{(k)}$is the $g$-th power sum of all the reciprocals ofthe $n-k$ zeros of $p(x)$ lying outside the disc$D(0,1/2)$. The leftmost equation of(\ref{equaA.4}) is verified by means of thedifferentiation of $p(x)=\prod_{j=1}^n(x-z_j)$.The middle equation of (\ref{equaA.4}) is impliedby the following decompositions, where $|x|=1$: $$\frac{1}{x-z_j}=\frac{1}{x}\sum_{h=0}^{\infty}(\frac{z_j}{x})^h\  for\ j\leq k,$$$$ \frac{1}{x-z_j}=-\frac{1}{z_j}\sum_{h=0}^{\infty}(\frac{x}{z_j})^h\for\ j> k.$$(\ref{equaA.3}) and (\ref{equaA.4}) together imply that$$ s_g^*=\sum_{i=-\infty}^{\infty} c_{-g-1+iQ},\ \ g=1,\ldots,Q-1.$$On the other hand, (\ref{equaA.4}) impliesthat $s_g^{(k)}=c_{-g-1}$ for $g \geq 0, \  s_g^{(k)}=-c_{-g-1}$ for $g<0$. We assumed in(\ref{equaA.3}) that $0<g<Q$. Therefore,$$s_g^*-s_g^{(k)} = \sum_{h=1}^{\infty}(s_{hQ+g}^{(k)}-s_{g-hQ}^{(k)}),$$\begin{equation}\label{equaA.5}|s_g^*-s_g^{(k)}|\leq \frac{kz^{Q+g}+(n-k)z^{Q-g}}{1-z^Q},\end{equation}where $0<g<Q$,$$z=\max_{1\leq j\leq n} \min \{|z_j|,\ 1/|z_j|\},$$and $z\leq 1/2$, since $i.r.(D)\geq 4$.\begin{remark}\label{remaA.1}The equations (\ref{equaA.3}) can be consideredan effective quadratureformulae for approximating the integrals$$ s_g^{(k)}=\frac{1}{2\pi\sqrt{-1}}\int_{\Gamma}(x^gp'(x)/p(x))dx$$ where$\Gamma$ denotes the unit circle $\{x: |x|=1\}$,$0<g<Q$, and (\ref{equaA.5})bounds the error of numerical integration.\end{remark} To compute the values $s_g^*$, we will performthree DFTs, each one on $Q$ points. Two DFTssuffice in order to compute $p(x)$ and $p'(x)$on the $Q$ nodes $\omega^q,\ q=0,1,\ldots, Q-1$,where $\omega=\exp(2\pi\sqrt{-1}/Q)$denotes a primitive $Q$-th root of 1. Thethird DFT computes the product of the Fouriermatrix [$\omega^{gh}$, $g=0,1,\ldots, Q-1$,$h=1,2,\ldots, n+1$] by the vector$p(\omega^q)/p'(\omega^q)$, $q=0,1,\ldots, Q-1$.$Q$ should be sufficiently large to guarantee thedesired precision of the approximation tothe power sums. Based on the error estimates(\ref{equaA.5}), Sch\"onhage in\cite{sc82} requires $Q$ to be of the order $n^2$,but this is only needed to handle the case where$i.r.(D)\geq1+1/O(n)$, which corresponds to$z=1-1/O(n)$ in (\ref{equaA.5}). In ourcase, $i.r.(D)\geq 4$; this corresponds to$z\leq 1/2$ in (\ref{equaA.5}). In this case,it suffices to choose $Q=O(n)$. Then the costof the approximate evaluation of all thepower sums $s_g^{(k)}$ is $O_A(n\log n)$. We are now prepared to deduce (\ref{equaA.1}) byfollowing \cite{sc82}, section~10 (compare \cite{sc82},section~11, or our Appendix B,for an alternative algorithm). At the cost$O_A(n\log n)$, we will obtain an initialapproximation $F_0(x)$ to the factor $F(x)$ of(\ref{equaA.2}) from approximations to thepower sums $s_g^{(k)}$ (see \cite{sc82}, section 13,or [P90], [BP94],solution of Problem~1.4.8) and then computean approximation $G_0(x)$ to the factor $G(x)$ of$p(x)$, by means of polynomial division, which can bereduced to inverting a triangular Toeplitz matrix(see [BP86] or [BP94]). The approximation errors are kept sufficiently small,so that the initial approximations to the factors can berapidly refined by means of Newton's iteration.Specifically, at Newton's iteration step $h$ we startwith the pair $F_{h-1}(x)$ and $G_{h-1}(x)$ andcompute the pair of polynomials $F_h(x)$ (monic), $G_h(x)$for $h=1,2,\ldots $, such that $\deg~F_h(x)=k$ monic,$\deg~G_h(x)=n-k$,\begin{equation}\label{equaA.6}\| F_h(x)G_h(x)-p(x)\|=\varepsilon_h,\end{equation}\begin{equation}\label{equaA.7}\varepsilon_{h+1}\leq \varepsilon_h^{1.5}.\end{equation}Here and hereafter, $\|\sum_i u_ix^i\|$ denotes$\sum_i|u_i|$, and with no lossof generality, we assume that $\|p(x)\|=1$. It isalso assumed that we have initially approximatedthe power sums $s_g^{(k)}$ and thencomputed the two polynomials $F_0(x)$ and$G_0(x)$ (see above), such that $F_0(x)$ is monic,\begin{equation}\label{equaA.8}\varepsilon_0\leq \eta/8;\ \  \varepsilon_0\leq(7\eta)^4/(k^2(7\eta+9k)^2 2^{2n+2k+2});\end{equation}\begin{equation}\label{equaA.9}\eta=\min_{|x|=1}|p(x)|\leq \|p(x)\|=1.\end{equation}Note that (\ref{equaA.7}) guarantees very rapidconvergence of $F_h(x)$ to $F(x)$ and $G_h(x)$ to$G(x)$, whereas (\ref{equaA.9}) and the secondbound of (\ref{equaA.8}) imply the first bound$\varepsilon_0\leq\eta/8$, since we may assume that$1\leq k<n$. We will start with recalling some knownauxiliary results (compare [Ah], [M] andsections~4 and 10 of \cite{sc82}). \begin{fact}\label{factA.1}For any polynomial $p(x)$, we have $$\|p(x)\|\geq \max_{|x|\leq 1}|p(x)|=\max_{|x|=1}|p(x)|.$$\end{fact} \begin{fact}\label{factA.2}$ 2^{d-1}\|q(x)r(x)\|\geq 2^{d-1}\max_{|x|=1}|q(x)r(x)|\geq\|q(x)\|\ \|r(x)\|$, for any pair of polynomials $q(x)$and $r(x)$ whose product has degree $d$.\end{fact} \begin{fact}\label{factA.3} Let$p^*(x)=h\prod_{i=1}^k(x-z_i^*)\\prod_{i=k+1}^n(1-x/z_i^*)$, where $|z_i^*|\leq 1$for $i\leq k;\ |z_i^*|\geq 1$for $i>k$. Then $|h|\geq \|p^*(x)\|/2^n$.\end{fact}In particular, for $p^*(x)=p(x)$, we obtain that$|h|>1/2^n$, since $\|p(x)\|=1$. Let us now deduce that \begin{equation}\label{equaA.9a}\eta\geq 1/2^{n+k},\end{equation}for $\eta$ of (\ref{equaA.9}),provided that $\|p(x)\|=1,\ |z_i|\leq 1/2$ for$i\leq k,\ |z_i|\geq 2$ for $i>k$. Due to thisbound on $\eta$, we may satisfy (\ref{equaA.8})by choosing $Q=Cn$ in (\ref{equaA.3}), for a sufficiently large constant $C$. To obtain (\ref{equaA.9a}), just observe that$$ |p(x)|=h\prod_{i=1}^k| x-z_i |\\prod_{i=k+1}^n|1-x/z_i |\geq2^{-n} (|x|-1/2)^k(1-|x|/2)^k=1/2^{2n}$$ provided that$|z_i|\leq 1/2$ for $i\leq k;\ |z_i|\geq 2$ for $i>k$, $|x|=1$.We will also use the following well-knowntheorem (see [Ah]):\begin{theorem}\label{theoA.1}{\rm (Rouch\'{e}).}If two functions $f_1(x)$ and $f_2(x)$ are analytic and bounded on a disc $D$ and if $|f_1(x)|>|f_2(x)|$on the boundary circle of the disc $D$, then bothfunctions $f_1(x)$ and $f_1(x)+f_2(x)$ have thesame number of zeros in the disc $D$.\end{theorem} Next, for a given pair of polynomials, $F_m(x)$ ofdegree $k$ and $G_m(x)$ of degree $n-k$, satisfying(\ref{equaA.6}) for $h=m$, we will compute twopolynomials, $F_{m+1}(x)$ of degree $k$ and$G_{m+1}(x)$ of degree $n-k$, such that$\varepsilon_{m+1}$ satisfying (\ref{equaA.6}) for$h=m+1$ also satisfies (\ref{equaA.7}) for$h=m$. To solve this task, we will followsection~10 of \cite{sc82}. We will firstcompute two polynomials $f_m(x)$ and $g_m(x)$ ofdegrees at most $k-1$ and at most $n-k-1$,respectively, such that\begin{equation}\label{equaA.11}\frac{p(x)-F_m(x)G_m(x)}{F_m(x)G_m(x)}=\frac{f_m(x)}{F_m(x)}+\frac{g_m(x)}{G_m(x)},\end{equation}and then we will compute\begin{equation}\label{equaA.12}F_{m+1}(x)=F_m(x)+f_m(x),\  \G_{m+1}(x)=G_m(x)+g_m(x).\end{equation}Since (\ref{equaA.11}) is an incomplete partialfraction decomposition, we may compute $f_m(x)$ and$g_m(x)$ by means of Euclid's extended algorithm,in $O(n\log n + l\log^2 l)$ operations,$l=\min\{k,\ n-k\}$ (see [AHU] or [BP94]). Our next task is to deduce (\ref{equaA.7}), for $h=m$.We assume that (\ref{equaA.7}) holds for$h=0,1,\ldots,m-1$, and since$\varepsilon_0\leq \eta/8\leq 1/8$ [see(\ref{equaA.8})], we have\begin{equation}\label{equaA.13}\varepsilon_m\leq \varepsilon_{m-1}\leq\varepsilon_0\leq\eta/8\leq 1/8.\end{equation}Now, from (\ref{equaA.12}), we obtain that\begin{eqnarray*}\lefteqn{F_{m+1}(x)G_{m+1}(x)-p(x)=} \\& & F_m(x)G_m(x)-p(x)+f_m(x)G_m(x)+F_m(x)g_m(x)+f_m(x)g_m(x).\end{eqnarray*}Combine this equation with (\ref{equaA.11}) to obtain that$$F_{m+1}(x)G_{m+1}(x)-p(x)=f_m(x)g_m(x)$$and, therefore,$$\varepsilon_{m+1}=\|F_{m+1}(x)G_{m+1}(x)-p(x)\|=\|f_m(x)g_m(x)\|\leq\|f_m(x)\|\ \|g_m(x)\|.$$We will also need the following auxiliary estimates from \cite{sc82},section~10 [where we use $\eta$of (\ref{equaA.9}) and $\varepsilon_h$ of (\ref{equaA.6}) for $h=m$]:\begin{lemma}\label{lemmA.1}$\|f_m(x)\|\leq (\varepsilon_m/(\eta-\varepsilon_m))\|F'_m(x)\|\leq(8/7)(\varepsilon_m/\eta)k\|F_m(x)\|$.\end{lemma}\begin{lemma}\label{lemmA.2}$\|g_m(x)\|\leq 2^{n+k-2}\varepsilon_m(1+(9/7)(k/\eta))/\|F_m(x)\|$.\end{lemma} {\em Proof of Lemma~\ref{lemmA.1}.} In the simple integral representation,$$ f_m(x)=\frac{1}{2\pi\sqrt{-1}}\int_{|t|=1}\frac{f_m(t)}{F_m(t)}\ \frac{F_m(x)-F_m(t)}{x-t}dt,$$first replace $\frac{f_m(t)}{F_m(t)}$ by $\frac{f_m(t)G_m(t)}{F_m(t)G_m(t)}$and then by $\frac{f_m(t)G_m(t)+g_m(t)F_m(t)}{F_m(t)G_m(t)}$. The latterreplacement is validated by the observation that$\int_{|t|=1}\frac{g_m(t)}{G_m(t)}\ \frac{F_m(x)-F_m(t)}{x-t}dt=0$, becausethe function $\frac{g_m(t)}{G_m(t)}\ \frac{F_m(x)-F_m(t)}{x-t}$is analytic in $t$, for $|t|\leq 1$ (see Lemma~\ref{lemmA.3}). Due to(\ref{equaA.11}), we have\begin{equation}\label{equaA.14}f_m(x)G_m(x)+g_m(x)F_m(x)=p(x)-F_m(x)G_m(x)\end{equation}and $$f_m(x)=\frac{1}{2\pi\sqrt{-1}}\int_{|t|=1}\frac{p(t)-F_m(t)G_m(t)}{F_m(t)G_m(t)}\ \frac{F_m(x)-F_m(t)}{x-t}dt.$$By applying this equation coefficient-wise, we obtain that $\|f_m(x)\|\leq\frac{\varepsilon_m}{\eta-\varepsilon_m}\|F'_m(x)\|$, which is the firstbound of Lemma~\ref{lemmA.1}. The second bound of Lemma~\ref{lemmA.1}follows since $\|F'_m(x)\|\leq k\|F_m(x)\|$ and $\varepsilon_m\leq \eta/8$[see (A.13)].  $\Box$ {\em Proof of Lemma~\ref{lemmA.2}.} Multiply both sides of (A.14) by$F_m(x)$ and obtain that $$g_m(x)F^2_m(x)=F_m(x)(p(x)-F_m(x)G_m(x))+f_m(x)(p(x)-F_m(x)G_m(x))-f_m(x)p(x),$$$$\|g_m(x)F^2_m(x)\|\leq \|F_m(x)\|\ \|p(x)-F_m(x)G_m(x)\|+\|f_m(x)\|\(\|p(x)-F_m(x)G_m(x)\|+\|p(x)\|).$$Now substitute the assumed equation $\|p(x)\|=1$, apply (\ref{equaA.6})for $h=m$, recall Lemma~\ref{lemmA.1}, and obtain that$$\|g_m(x)F^2_m(x)\|\leq \|F_m(x)\|\varepsilon_m +(8/7)(\varepsilon_m/\eta)k\|F_m(x)\|(\varepsilon_m+1).$$By combining the bound of Fact~\ref{factA.2} for $q(x)=g_m(x)$,$r(x)=F_m^2(x)$, $d=n+k-1$, with the latter upper bound on$\|g_m(x)F_m^2(x)\|$ and with (A.13), we arrive atLemma~\ref{lemmA.2}.\ \  $\Box$ To complete the proof of Lemma~\ref{lemmA.1} we need\begin{lemma}\label{lemmA.3} For all $h,\ h=0,1,\ldots$, allthe zeros of $F_h(x)$ lie strictly inside the unit disc$D(0,1)$, and all the zeros of $G_h(x)$ lie outside this disc.\end{lemma} {\em Proof.} The second assertion [about $G_h(x)$] followsfrom the first one [about $F_h(x)$] and from (Rouch\'{e}'s)Theorem~\ref{theoA.1};indeed, due to this theorem and to the bounds(\ref{equaA.6}) and (\ref{equaA.8}), the polynomials$p(x)$ and $F_h(x)G_h(x)$ have the same numberof zeros in $D(0,1)$. The first assertion is initiallyensured for $F_0(x)$ [by the definition of $F_0(x)$] andthen is recursively extended from $F_h(x)$ to $F_{h+1}(x)$,for $h=0,1,\ldots$, by applying (Rouch\'{e}'s)Theorem~\ref{theoA.1} to $F_h(x)G_h(x)$ and$F_{h+1}(x)G_h(x)=F_h(x)G_h(x)+f_h(x)G_h(x)$.We only need to prove that$\gamma_h=|f_h(x)G_h(x)|<|F_h(x)G_h(x)|$, for all$x$ satisfying $|x|=1$ and for $h=0,1,\ldots$, in order to justifythe application of (Rouch\'{e}'s) Theorem~\ref{theoA.1}.We have $$\gamma_h \leq\|f_h(x)\|\  \|G_h(x)\|.$$Apply Lemma~\ref{lemmA.1} for $m=h$ (this is valid, sincewe are now aiming at $F_{h+1}(x)$, in our inductiveextension of Lemma~\ref{lemmA.3} from $h$to $h+1$) and deduce that $$\gamma_h\leq (8/7)(\varepsilon_h/\eta)k\|F_h(x)\|\ \|G_h(x)\|.$$Apply Fact~\ref{factA.2} and obtain that$$ \gamma_h\leq(8/7)(\varepsilon_h/\eta)k 2^{n-1}\|F_h(x)G_h(x)\|.$$by using (\ref{equaA.6}) and (A.13) for $m=h$, we nowdeduce that $$\|F_h(x)G_h(x)\|\leq \|p(x)\|+\varepsilon_h\leq 1+\varepsilon_h\leq 9/8,$$ and therefore, $$\gamma_h\leq (9/7)(\varepsilon_0/\eta)k2^{n-1}\leq(9/7)7^4\eta^3/(k(7\eta +9k)^22^{n+3+2k})<$$$$7\eta /8=\eta-\varepsilon_h\leq |p(x)|-\|p(x)-F_h(x)G_h(x)\|.$$Apply Fact~\ref{factA.1} and deduce that$$\gamma_h<|p(x)|-|p(x)-F_h(x)G_h(x)|\leq |F_h(x)G_h(x)|,$$for all $x$ satisfying $|x|=1$.\ \  $\Box$ If we complemented (\ref{equaA.7}) and (\ref{equaA.8}) by the requirement that$\varepsilon_h<2^{-7n},\  \ h=0,1,\ldots,$ then we could have skippedthe above proof and immediately deduced Lemma~\ref{lemmA.3}from the following extension of Ostrowski's theorem:\begin{theorem}\label{theoA.2}[Sc85].Let $$p(x)=p_n\prod_{j=1}^n (x-z_j),\ \ p^*(x)=p^*_n\prod_{j=1}^n (x-z_j^*),$$$$\|p^*(x)-p(x)\|\leq q\|p(x)\|,\ \ q<2^{-7n},$$$$|z_j|\leq 1,\ \ j=1,\ldots,k;\ \  |z_j|\geq 1,\ \ j=k+1,\ldots,n.$$ Then$$ |z_j^*-z_j|<9\sqrt[n]q,\ \ j=1,\ldots,k; $$$$ |1/z_j^*-1/z_j|<9\sqrt[n]q,\ \ j=k+1,\ldots,n. $$\end{theorem} Next, by combining Lemmas~\ref{lemmA.1} and \ref{lemmA.2},we obtain that$$\varepsilon_{m+1} = \|f_m(x)\|\ \|g_m(x)\|\leq \varepsilon_m^2k2^{n+k-2}\frac{8}{7\eta}(1+\frac{9k}{7\eta}),$$ and therefore, $$\varepsilon_{m+1}\leq \varepsilon_m^22^{n+k+1}\frac{(7\eta + 9k) k}{49\eta^2}.$$By combining this bound with (\ref{equaA.8}), we obtain(\ref{equaA.7}), for $h=m$. Due to (\ref{equaA.7}), we ensure splitting with the accuracy$\varepsilon_h$ in $O(\log \log (\varepsilon_0/\varepsilon_h))$ steps,where each step is performed by using$O(n\log n + l\log^2 l)$ operations, for$l=\min\{k,\ n-k\}$ [since the complexity of eachstep is dominated by the complexity ofcomputing the partial fraction decomposition(\ref{equaA.11})]. Note that these$O(\log \log (\varepsilon_0/\varepsilon_h))$ steps of theinner loop are included into the outer loop of $n-1$ stepsof recursive splitting. If we seek splitting of $p(x)$ intolinear factors within the error bound $E(p)$, then to achieveour goal, we should require that in $u-1$ steps of the outer loopof the recursive splitting, the $u$ computed factors$P_1(x), P_2(x), \ldots, P_u(x)$ satisfy the bounds\begin{equation}\label{equaA.15}E_u=\|p(x)-\prod_{i=1}^uP_i(x)\|\leq uE(p)\|p(x)\|/n,\end{equation}for $u=2,\ldots,n$(see \cite{sc82}).Let us now examine what this meansin terms of splitting a factor $P_j(x)$ into two factors$\hat{F}(x)$ and $\hat{G}(x)$ at the $l$-th step for$1\leq l<n$, $1\leq j\leq l$. Let, say, $j=l$ and$$ \|P_l(x)-\hat{F}(x)\hat{G}(x)\|=e_l\|P_l(x)\|.$$Combining the latter equation with (\ref{equaA.15})for $u=l$ gives us that$$ \begin{array}{lcl}E_{l+1}&=&\|p(x)-\hat{F}(x)\hat{G}(x)\prod_{i=1}^{l-1}P_i(x)\|\\&\leq & E_l + \|\prod_{i=1}^lP_i(x)-\hat{F}(x)\hat{G}(x)\prod_{i=1}^{l-1}P_i(x)\|\\&\leq & lE(p)\|p(x)\|/n +e_l\|P_l(x)\|\\|\prod_{i=1}^{l-1} P_i(x)\|\\ &\leq & lE(p) \|p(x)\|/n+ e_l\prod_{i=1}^l \|P_i(x)\|. \end{array} $$Now we can see that it is sufficient to impose the bound\begin{equation}\label{equaA.16}e_l\leq E(p)\|p(x)\|/(n\prod_{i=1}^l\|P_i(x)\|)\end{equation}in order to arrive at (\ref{equaA.15}) for $u=l+1$. By first applying Fact~\ref{factA.2} and then(\ref{equaA.15}) for $u=l$, we obtain that $$\prod_{i=1}^l\|P_i(x)\|\leq 2^{n-1}\|\prod_{i=1}^lP_i(x)\|\leq 2^{n-1}(1+lE(p)/n)\|p(x)\|,$$which does not exceed $2^n\|p(x)\|$ for $l\leq n$,$E(p)\leq 1$. Therefore, it suffices to require that\begin{equation}\label{equaA.17}e_l\leq E(p)/(n2^n)\ \  {\rm for} \ \ {\rm all} \ \ l \ ,\end{equation}in order to ensure (\ref{equaA.16}) and therefore(\ref{equaA.15}) for all $u$. Let us now estimate the computational complexity of theentire algorithm under the latter requirement.We recall the estimates of Ostrowski's theorem,according to which (\ref{equaA.15}) for$u=n$ implies approximating the zeros of$p(x)$ by the zeros of $F_h(x)$ and $G_h(x)$within the errors $O((E(p))^{1/n})$(for some refined error estimates, applyTheorem~\ref{theoA.2} if$E(p)<2^{-7n}$ or see \cite{sc82},section~19 ).  By using (\ref{equaA.17}) wetranslate the bounds on $E(p)$ into ones on $e_l$ and$\varepsilon_h$ for all $l$ and $h$ and obtainthat the choice of $\varepsilon_h$ satisfying$\log \log (\varepsilon_0/\varepsilon_h)=O(\log(bn))$suffices for our purpose. Now recall that $\varepsilon_h\leq \varepsilon$, so thatthe outer loop requires $O(\log(bn))$ callsfor the inner loop. Since the complexity of each applicationof the inner loop is $O(l\log^2 l)$, it followsthat the overall complexityof the outer loop is bounded by (\ref{equaA.1})if we only consider the stages, where $l\leq n/\log n$.Each other stage splits a polynomial of degree $u$into $l$ nonconstant factors, each of degree at most$u-l+1\leq u-n/\log n$. Therefore, there can be at most$\log n$ recursive splittings with $l>n/\log n$, andsince in all cases $l\leq n$, their total contributionto the overall complexity bound is$O(n\log^3 n\log (bn))$, which is dominated bythe desired estimate (\ref{equaA.1}). \begin{remark}\label{remaA.2} The bound (\ref{equaA.1}) isachieved for all input instances $p(x)$. For a large classof input polynomials, all splittingsgive factors $F(x)$ and $G(x)$ of roughly the same orcomparable degrees, and in this case the cost ofthe computations decreases by roughly the factor of$n$ [P94], similarly to the phenomenon known,for instance, in the QUICKSORT algorithm.Our results of this paper facilitate achieving the balanced splitting byshifting to a neighboring polynomial $p(x) +x^i \Delta$.\end{remark}\begin{remark}\label{remaA.3} The computation of thepartial fraction decomposition (\ref{equaA.11}) requiresa higher precision, which suggests that the competitiveapproaches (of Appendix B)may still be superior, in terms of their Boolean (bit)complexity. While applying the algorithm ofthis section, however, one may use the various techniquesof symbolic computing [such as the Chinese remainderalgorithm and $p$-adic (Newton-Hensel's) lifting] todecrease the precision of computing at the stage ofcomputing the decomposition (\ref{equaA.11}).\end{remark}\vspace{.3in} \section{ \normalsize\bfAn Alternative Algorithm for Recursive Splitting intoFactors}\label{appB}\setcounter{equation}{0} Here we will recall (and slightly simplify by removing somecomplicated contour integration techniques) an alternativealgorithm, of section~11 of \cite{sc82}, which evaluates thepolynomials $f_m(x)$ and $g_m(x)$ satisfying(A.11) and (A.12), for two given polynomials$F_m(x)$ and $G_m(x)$, approximating two factors of $p(x)=F(x)G(x)$.In its original form, this alternativealgorithm only supports the arithmetic complexity estimate of$O(n^3\log n \log (bn))$; our modificationsupports the estimate (\ref{equaA.1}) and Proposition~\ref{propA.1};furthermore, according to a claim in \cite{sc82}, thisalgorithm leads to decreasing (by roughly the factor $\log n$)the upper bound on the number of bit-operations involved inthis stage (compare Remark~\ref{remaA.3} inour Appendix A). To describe the algorithm of \cite{sc82},with our modification, we first assume thatwe are given some auxiliary polynomial $H_m(x)$ of degreeless than $k$ and such that\begin{equation}\label{equaf.1}H_m(x)G_m(x)+J_m(x)F_m(x)=1\end{equation}for some polynomial $J_m(x)$. Therefore,$$ H_m(x)(p(x)-F_m(x)G_m(x))G_m(x)+J_m(x)(p(x)-F_m(x)G_m(x))F_m(x)=p(x)-F_m(x)G_m(x).$$Comparison with (\ref{equaA.12}) suggests the choice of$(p(x)-F_m(x)G_m(x))H_m(x)$ as $f_m(x)$ and of$(p(x)-F_m(x)G_m(x))J_m(x)$ as $g_m(x)$, provided  that we bound the degrees of $f_m(x)$ by $k-1$and of $g_m(x)$ by $n-k-1$. Thus, having the desiredpolynomial $H_m(x)$ of (\ref{equaf.1}), one should first compute\begin{equation}\label{equaf.2}f_m(x)=p(x)H_m(x) \bmod F_m(x),\end{equation}so that $\deg f_m(x)< \deg F_m(x)=k$, and then obtain $g_m(x)$as the quotient of the division of the polynomial$p(x)-F_m(x)G_m(x)-f_m(x)G_m(x)$ by $F_m(x)$, so that$\deg g_m(x)=\deg (p(x)-F_m(x)G_m(x)-f_m(x)G_m(x))-\deg F_m(x)<n-k$. Now suppose that, instead of a polynomial $H_m(x)$ satisfying(\ref{equaf.1}), we have an approximation $H_m^*(x)$ to $H_m(x)$ such that\begin{equation}\label{equaf.3}H_m^*(x)G_m(x) + J_m(x)F_m(x)=1-D_m(x),\ \ \deg H_m^*(x)<k,\\deg D_m(x)<n,\end{equation}\begin{equation}\label{equaf.4}\|p(x)-F_m(x)G_m(x)\|=\varepsilon_m\leq \delta_m; \ \ \|D_m(x)\|\leq \delta_m;\end{equation}\begin{equation}\label{equaf.5}\delta_m=\delta_0^{1.5^{m}},\  \\delta_0\leq \eta^4/(k^42^{3n+k+1})\leq \eta/256,\end{equation}for $\eta$ of (\ref{equaA.9}), $\eta\leq 1$, where in(\ref{equaf.5})  the second inequality follows from thefirst one, since $\eta\leq 1\leq k<n$.Then the algorithm of \cite{sc82} first computes\begin{equation}\label{equaf.6}f_m^*(x)=p(x)H_m^*(x) \bmod F_m(x) =(p(x)-F_m(x)G_m(x))H^*_m(x) \bmod F_m(x),\end{equation}via polynomial multiplication modulo $F_m(x)$ [compare(\ref{equaf.2})], then$g_m^*(x)$ (by means of polynomial division) such that\begin{equation}\label{equaf.7}p(x)-F_m(x)G_m(x) -f_m^*(x)G_m(x)=g_m^*(x)F_m(x)+r_m(x),\end{equation}$\deg r_m(x)<k$ [where we assume that$\deg r_m(x)=-\infty$ if $r_m(x)=0$], and finally $$F_{m+1}^*(x)=F_m(x)+f^*_m(x),\ \ G_{m+1}^*(x)=G_m(x)+g^*_m(x).$$Given $H_m^*(x)$, this computation of $F_{m+1}^*(x)$and $G_{m+1}^*(x)$only involves $O (n\log n)$ arithmetic operations. Our next goal is to estimate that\begin{equation}\label{equaf.8}\varepsilon^*_{m+1}=\|p(x)-F_{m+1}^*(x)G_{m+1}^*(x)\|\leq \delta_{m+1}=\delta_m^{1.5}\end{equation}provided that (\ref{equaf.3})--(\ref{equaf.7}) hold.It follows from (\ref{equaf.7}) that\begin{eqnarray}\varepsilon^*_{m+1}&=&\|p(x)-F_{m+1}^*(x)G_{m+1}^*(x)\|\nonumber \\  &\leq &\|f_m^*(x)g_m^*(x)-r(x)\|\nonumber \\   &\leq &\|f_m^*(x)\|\ \|g_m^*(x)\|+\|r_m(x)\|.\label{equaf.9}\end{eqnarray}Due to (\ref{equaf.3}) and (\ref{equaf.6}), we have $$f_m^*(x)G_m(x)=(p(x)-F_m(x)G_m(x))(1-D_m(x))\bmod F_m(x).$$In our further analysis, we will use the following extension ofLemma~\ref{lemmA.1} (whose proof is similar to one ofLemma~\ref{lemmA.1}):\begin{lemma}\label{lemmf.1}(\cite{sc82}, section~10). Let$s(x)=u(x)g(x)+v(x)f(x)$ for some polynomials$f(x),\  g(x),\  s(x),\  u(x)$ and $v(x)$ such that$\deg u(x)< \deg f(x)$. Let $f(x)$ have exactly $k$zeros, all lying in $D(0,1/2)$, and let $g(x)$ have exactly $n-k$zeros, all lying outside $D(0,2)$. Let $|f(x)g(x)|\geq \eta^* $ if$|x|=1$. Then$$ \|u(x)\|\leq \|s(x)\|\  \\|f'(x)\|/\eta^*\leq k\|s(x)\|\ \|f(x)\|/\eta^*.$$\end{lemma}We will several times apply Lemma~\ref{lemmf.1} for$f(x)=F_m(x),\ g(x)=G_m(x),\ \eta^*=\eta-\delta_m$[see (\ref{equaA.9}) and (\ref{equaf.5})], so that$\eta^*\geq 255\eta/256$. In the first application,we will also assume that$u(x)=f_m^*(x),\ s(x)=(p(x)-F_m(x)G_m(x))(1-D_m(x))$.Taking into account (\ref{equaf.4})--(\ref{equaf.6}),we now deduce that\begin{equation}\label{equaf.10}\|f_m^*(x)\|\leq (256/255)(k/\eta)\delta_m (1+\delta_m)\|F_m(x)\|.\end{equation}Similarly, (\ref{equaf.3}), (\ref{equaf.6}) and (\ref{equaf.7})together imply that $$r_m(x)G_m(x)=D_m(x)(p(x)-F_m(x)G_m(x))G_m(x)\bmod F_m(x). $$We again apply Lemma~\ref{lemmf.1}, for$f(x)=F_m(x),\ g(x)=G_m(x),\ \eta^*=\eta-\delta_m$, but thistime set $u(x)=r_m(x),\ s(x)=D_m(x)(p(x)-F_m(x)G_m(x))G_m(x)$and deduce that\begin{eqnarray}\|r_m(x)\| & \leq & (256/255)(k/\eta)\delta_m^2 \|F_m(x)\|\\|G_m(x)\|\nonumber\\  & \leq & (256/255)(k/\eta)\delta_m^2  2^{n-1}(1+\delta_m).\label{equaf.11}\end{eqnarray}Here again we take into account (\ref{equaf.4})and (\ref{equaf.5}) in order to deduce theformer inequality of (\ref{equaf.11}), whose latterinequality then follows from Fact~\ref{factA.2} for$q(x)=F_m(x),\  r(x)=G_m(x)$ and from the simplebound $\|F_m(x)G_m(x)\|\leq\|p(x)\|+\delta_m=1+\delta_m$,implied by (\ref{equaf.4}).From (\ref{equaf.7}) one also deduces that $$\|g_m^*(x)F_m^2(x)\|\leq \|r_m(x)\|\ \|F_m(x)\| +\delta_m\|F_m(x)\| +\|f_m^*(x)\|\ \|F_m(x)G_m(x)\|.$$Then, by applying Fact~\ref{factA.2} again, oneobtains that$$ \|g_m^*(x)\|\leq 2^{n+k-2}\|g_m^*(x)F_m^2(x)\|/\|F_m(x)\|^2\leq $$\[2^{n+k-2}\left(\|r_m(x)\|+\delta_m+\|f_m^*(x)\|(1+\delta_m)/\|F_m(x)\|\right)/\|F_m(x)\|.\]We now combine the above estimate with (\ref{equaf.10})and (\ref{equaf.11})  and deduce that$$ \|g_m^*(x)\|\leq [1+(256 k/(255\eta ))(1+\delta_m)(\delta_m2^{n-1}+(1+\delta_m))]2^{n+k-2}\delta_m/\|F_m(x)\|.$$Since $\eta\leq 1\leq k<n$, we easily deduce from the latterbound and from (\ref{equaf.5}) that\begin{equation}\label{equaf.12}\|g_m^*(x)\|\leq (0.51)2^{n+k}\delta_mk/(\eta\|F_m(x)\|).\end{equation}By combining (\ref{equaf.9})--(\ref{equaf.12}), we nowobtain that $$\begin{array}{lll}\varepsilon^*_{m+1}& = & \|p(x)-F_{m+1}^*(x)G_{m+1}^*(x)\|\\  & \leq & (256/255)\delta_m^2(1+\delta_m)(k/\eta)(((0.51)2^{n+k}k/\eta)+2^{n-1})\\  &\leq & 2^{n+k}\delta_m^2(k/\eta)^2\end{array} $$since $(256/255)(1+\delta_m)(0.51+\eta/(k2^{k+1}))<1$.From the latter boundon $\varepsilon^*_{m+1}$, we immediately deduce (\ref{equaf.8}). It remains to ensure (\ref{equaf.8}) for $m=0$ and for$\delta_0$ that satisfies(\ref{equaf.5}) and then recursively apply (\ref{equaf.8})for $m=1,2,\ldots$. In this way, we rapidly improve theinitial approximations to $F(x)$ by $F_m(x)$  andto $G(x)$ by $G_m(x)$, and it follows that we mayapproximate the splitting of$p(x)$ into the factors $F(x)$ and $G(x)$at the cost bounded according to Proposition~\ref{propA.1} if we ignore the cost of obtaining a polynomial $H_m^*(x)$ thatsatisfies the relations (\ref{equaf.3})--(\ref{equaf.5}). Thus, it remains to specify the algorithms and complexityfor the following operations:\begin{description}\item[a)] the initial choice of $H_0^*(x)$ [which will also  give us $D_0(x)=1-H_0^*(x)G_0(x) \bmod F_0(x)$, due to (\ref{equaf.3}) applied for $m=0$] and \item[b)] the transition from $H_m^*(x)$ to $H_{m+1}^*(x)$,for $m=0,1,\ldots$.\end{description}We will start with the harder part a). The polynomial$H_0^*(x)$ will be computed as a numerical approximation tothe polynomial $H_0(x)$defined by the partial fraction decomposition$$ \frac{1}{F_0(x)G_0(x)}=\frac{H_0(x)}{F_0(x)}+\frac{J_0(x)}{G_0(x)},\ \\deg J_0(x)<n-k$$[compare (\ref{equaf.1})] andhaving the following integral representation(compare \cite{sc82}, section~12):$$H_0(x)=\frac{1}{2\pi\sqrt{-1}}\int_{|t|=1}\frac{F_0(x)-F_0(t)}{F_0(t)G_0(t)(x-t)}dt=\sum_{h=0}^{k-1}x^h\sum_{i=1+h}^kf_iv_{i-h-1},$$where $$F_0(x)=\sum_{i=0}^kf_ix^i,\ \v_h=\frac{1}{2\pi\sqrt{-1}}\int_{|t|=1}\frac{t^hdt}{F_0(t)G_0(t)},\ \h=0,1,\ldots,k.$$The values $v_h$ can be well approximated by using DFTs[compare (\ref{equaA.3})--(\ref{equaA.5})], since\begin{equation}\label{equaf.13} v_h=c_{-h-1}^*,\ \v_{p-1}^*(Q)=\frac{1}{Q}\sum_{j=0}^{Q-1}\frac{1}{F_0(\omega^j)G_0(\omega^j)}\omega^{pj}=\sum_{i=-\infty}^{\infty}c_{-p+iQ}^*,\end{equation}where \begin{equation}\label{equaf.14} \frac{1}{F_0(x)G_0(x)}=\sum_{i=-\infty}^{\infty}c_i^*x^i\end{equation}and where $v_h^*(Q)$ denote the computedapproximations to $v_h$ for all $h$. Then theconvolution gives us the desired coefficients$\sum_{i=1+h}^kf_iv_{i-h-1}^*(Q)$ of the polynomial$H_0^*(x)$ approximating $H_0(x)$.Clearly, the FFT-based complexity of computing$v_h^*=v_h(Q)$ for $h=0,\ldots,k$ is bounded by$O_A(\log Q, n+Q)$. Next, we will prove that the desired bound$\|D_0(x)\|\leq \delta_0$, for $\delta_0$ of(\ref{equaf.5}), can be achieved already where$Q=O(n)$, which will imply the overallcomplexity $O_A(\log n, n)$ for the computations ofpart a), that is, for approximating thepolynomial $H_0(x)$. We begin this proof by deducingfrom Fact~\ref{factA.2} that$\|F_0(x)\|\leq 2^{n-1}\|F_0(x)G_0(x)\|/\|G_0(x)\|$. Now,since $\|F_0(x)G_0(x)\|\leq\|p(x)\|+\delta_0\leq 1+1/256<2$,we have $ \|F_0(x)\|\leq 2^n/\|G_0(x)\|$, and we obtain that$$\|H_0^*(x)-H_0(x)\|\leq k\|F_0(x)\|\Delta(v)\leqk2^n\Delta(v)/\|G_0(x)\|,$$where $\Delta(v)=\Delta_Q(v)=\max_{0\leq h \leq k}|v_h-v_h^*(Q)|$. By comparing (\ref{equaf.1}) with (\ref{equaf.3}),we deduce that the bound $\|D_0(x)\|\leq \delta_0$ will beensured for $\delta_0$ satisfying (\ref{equaf.5}) if\begin{equation}\label{equaf.15}\Delta(v)\leq \eta^4/(k^52^{4n+k+1}).\end{equation}Thus, it remains to show that (\ref{equaf.15}) can be satisfied for$Q=O(n)$. From (\ref{equaf.13}) and (\ref{equaf.14}) we obtain that\begin{equation}\label{equaf.16}\Delta(v)\leq \max_{0\leq h\leq k}|(\sum_{i=-\infty}^{\infty}c_{-h+iQ}^*)-c_{-h}^*|.\end{equation}Thus, to solve our task, we only need to estimate$|c_s^*|$ from above, for $s\geq Q-k$ and $s\leq -Q-k$ .We will rely on\begin{proposition}\label{propf.1}(\cite{sc82}, sect.~9).Let $p(x),\ p^*(x)$ be two polynomialsof degrees at most $n$, $|p(x)|\geq \eta$ for$|x|=1,\ p(x)\not = 0$ for $1/q^2\leq |x|\leq q^2,\ q>1$.Then $$|p^*(x)|\geq \eta/(q+1)^n-\|p^*(x)-p(x)\|,$$ for$1/q\leq |x|\leq q$.\end{proposition}By applying Proposition~\ref{propf.1} for$p^*(x)=F_0(x)G_0(x)$, we deduce that\begin{equation}\label{equaf.17}|F_0(x)G_0(x)|\geq \eta/(q+1)^n-\delta_0 \ \ {\rm for}\ \1/q\leq |x|\leq q.\end{equation}We may set $q=\sqrt 2$ since, by assumption, the unitdisc $D(0,1/2)$ is 4-isolated.We now recall that $1\leq k<n$,$1/2^{n+k}\leq \eta\leq 1$ [see (\ref{equaA.9a})],also recall (\ref{equaf.4}) and (\ref{equaf.5}) anddeduce from (\ref{equaf.17}) that$|F_0(x)G_0(x)|\geq (3^{-n}-2^{-3n})\eta >3^{-n}\eta/2 =\eta 2^{-n\log 3 -1}> 2^{-2n-1}\eta \geq 2^{-k-3n-1}$, and therefore,\begin{equation}\label{equaf.18}1/|F_0(x)G_0(x)|<2^{k+3n+1},\ \  {\rm for}\ \  1/\sqrt 2\leq |x|\leq \sqrt 2.\end{equation}Now we divide both sides of (\ref{equaf.14}) by $x^{i+1}$,apply a well-known result from the analytic function theory,[Ah], and obtain that $$c_s^*=\frac{1}{2\pi\sqrt{-1}}\int \frac {x^{-s-1}\ dx}{F_0(x)G_0(x)},$$where the integration is along the boundaries of thetwo circles, $\{ x:|x|=1\}$ and $\{ x:|x|=2^{0.5\  {\rm sign}(s)}\}$, sign$(s) = 1$ if $s>0$, sign$(0) =0$, sign$(s)=-1$ if $s<0$. Dueto (\ref{equaf.18}), it follows that $|c_s^*|\leq 2^{-0.5 |s|+k+3n+2}$, for $s=0, \pm1,\pm 2,\ldots$.We may now rewrite (\ref{equaf.16}) as follows:$$ \Delta (v)\leq \sum_{i=1}^{\infty}|c_{-h+iQ}^*|+\sum_{i=-1}^{-\infty}|c_{-h+iQ}^*| \leq \sum_{i=1}^{\infty}(2^{0.5(h-iQ)+k+3n+2} +2^{-0.5(h+iQ)+k+3n+2}),$$ for $h=0,1,\ldots, k$. Therefore,for $Q>1$, we have:$$ \Delta(v)< 2^{3-0.5Q+1.5k+3n}.$$ We now apply(\ref{equaA.9a}) and arrive at(\ref{equaf.15}) already for$Q=22n+13k+10\log k +8$. This completes thecomplexity analysis of part  a), that is,of the stage of approximating $H_0(x)$. Let us now consider part b), that is, the transitionfrom $H_m^*(x)$ to $H_{m+1}^*(x)$, for $m=0,1,\ldots$,where we need to maintain the relations(\ref{equaf.3})--(\ref{equaf.5}) for all $m$. We will start with anauxilary result.\begin{lemma}\label{lemmf.2} Let polynomials$H_m^*(x)$, $D_m(x)$, $f_m^*(x)$, $g_m^*(x)$, $G_m(x)$,$F_m(x)$, $J_m(x)$ satisfy the relations(\ref{equaf.3})--(\ref{equaf.7}). Let\begin{equation}\label{equaf.21}D_{m,0}(x)=D_m(x)-H_m^*(x)g_m^*(x)-J_m(x)f_m^*(x).\end{equation}Then\begin{equation}\label{equaf.22}\deg D_{m,0}(x)<n,\end{equation}\begin{equation}\label{equaf.23}H_m^*(x)G_{m+1}^*(x)+J_m(x)F_{m+1}^*(x)=1-D_{m,0}(x),\end{equation}where $F_{m+1}^*(x)=F_m(x)+f_m^*(x)$,$G_{m+1}^*(x)=G_m(x)+g_m^*(x)$, and\begin{equation}\label{equaf.24}\|D_{m,0}(x)\|<2.1\delta_m(k/\eta)^2 2^{n+k}.\end{equation}\end{lemma}{\em Proof.} (\ref{equaf.3}) and(\ref{equaf.21}) together immediately imply(\ref{equaf.22}) and (\ref{equaf.23}). To deduce(\ref{equaf.24}), first apply Lemma~\ref{lemmf.1} to(\ref{equaf.3}) for $s(x)=1-D_m(x)$, $u(x)=H_m^*(x)$,$v(x)=J_m(x)$, $f(x)=F_m(x)$, and $g(x)=G_m(x)$.This application gives us the bound\begin{equation}\label{equaf.25}\|H_m^*(x)\|/\|F_m(x)\|\leq \|1-D_m(x)\|(k/\eta^*)\leq(1+\delta_m)256k/(255\eta).\end{equation}Then multiply the equation of (\ref{equaf.3}) by $F_m(x)$ anddeduce that $$\|J_m(x)F^2_m(x)\|\leq (1+\|D_m(x)\|)\|F_m(x)\|+\|H_m^*(x)\|\ \|F_m(x)G_m(x)\|.$$Substitute the bounds of (\ref{equaf.4}) and(\ref{equaf.25}), recall the assumptionthat $\|p(x)\|=1$ and obtain that $$\|J_m(x)F_m^2(x)\|\leq (1+\delta_m)\|F_m(x)\|(1+(1+\delta_m)256k/(255\eta))<2.01(1+\delta_m)^2\|F_m(x)\|k/\eta.$$Apply Fact~\ref{factA.2} and obtain that $$\|J_m(x)\|\ \|F_m^2(x)\|\leq 2^{n+k-2}\|J_m(x)F_m^2(x)\|$$and, consequently,\begin{equation}\label{equaf.26}\|J_m(x)\|\ \|F_m(x)\|\leq (2.01/4)2^{n+k}(1+\delta_m)^2k/\eta.\end{equation}Now combine the relations  (\ref{equaf.4}),(\ref{equaf.10}), (\ref{equaf.12}), (\ref{equaf.21}),(\ref{equaf.25}) and (\ref{equaf.26}) and obtain that\begin{eqnarray}\label{equaf.27}\|D_{m,0}(x)\|&\leq & \delta_m+(0.51)2^{n+k}\delta_m(k/\eta)^2(1+\delta_m)(256/255)+ \nonumber  \\ &  &(256/255)(k/\eta)^2\delta_m(1+\delta_m)^3(2.01/4)2^{n+k}\nonumber\\ &=&\delta_m(1+2^{n+k}(1+\delta_m)(k/\eta)^2(256/255)(0.51+(1+\delta_m)^22.01/4)) \nonumber \\&<& 2.1\delta_m2^{n+k}(k/\eta)^2,\nonumber\end{eqnarray}which completes the proof ofLemma~\ref{lemmf.2}.  $\Box$ Due to Lemma~\ref{lemmf.2} our task of part b) is narrowed:we only need to make transition from the pair of polynomials$H_m^*(x)$, $D_{m,0}(x)$ satisfying(\ref{equaf.21})--(\ref{equaf.24}) to the pair$H_{m+1}^*(x)$, $D_{m+1}(x)$, satisfying(\ref{equaf.3})--(\ref{equaf.5}) for $m$ replaced by$m+1$. To achieve this, we set $H_{m,0}(x)=H_m^*(x)$and then compute a sequence$\{ H_{m,i+1}(x), D_{m,i+1}(x) \}, \ \ i=0,1,\ldots$, ofthe auxiliary polynomials defined by the followingformulae:\begin{eqnarray}H_{m,i+1}(x)&=&H_{m,i}(x)(1+D_{m,i}(x))\bmod F_{m+1}^*(x),\label{equaf.28a}\\D_{m,i+1}(x)&=&D^2_{m,i}(x)\bmod F_{m+1}^*(x),\label{equaf.28}\end{eqnarray}$$\deg H_{m,i+1}(x)<k,\ \ i=0,1,\ldots.$$ We stop the computation and set$H_{m+1}^*(x)=H_{m,i}(x)$, $D_{m+1}(x)=D_{m,i}(x)$ as soon aswe satisfy (\ref{equaf.3})--(\ref{equaf.5})  forsome $i$ and for $m$ replaced by $m+1$. We need to extend the properties(\ref{equaf.22})--(\ref{equaf.24}) to $H_{m,i}(x)$ and$D_{m,i}(x)$. Let us now inductively extend (\ref{equaf.23})by deducing that\begin{equation}\label{equaf.27n}H_{m,j}(x)G_{m+1}^*(x)=1-D_{m,j}(x)\bmod F_{m+1}^*(x)\ \ {\rm for }\ \j=0,1,\ldots.\end{equation}Indeed, (\ref{equaf.23}) gives us this equation for $j=0$.For the transition from $j=i$ to $j=i+1$, multiplythe equation of (\ref{equaf.28a}) by $G_{m+1}^*(x)$ and thensubstitute (\ref{equaf.27n}) for $j=i$ on the right hand sideto obtain that $$\begin{array}{lcl}H_{m,i+1}(x)G_{m+1}^*(x)&=&H_{m,i}(x)G_{m+1}^*(x)(1+D_{m,i}(x))\bmod F_{m+1}^*(x)\\ &=&1-D_{m,i}^2(x)\bmod F_{m+1}^*(x)\ . \end{array} $$Now substitute the equation of (\ref{equaf.28}) on theright-hand side and arrive at (\ref{equaf.27n}) for $j=i+1$,which completes the inductive proof of (\ref{equaf.27n}). Next, we will extend the bound (\ref{equaf.24}).From (\ref{equaf.28}), we deduce that $$D_{m,i+1}(x)G_{m+1}^*(x)+J_{m,i}(x)F_{m+1}^*(x)=D_{m,i}^2(x)G_{m+1}^*(x) \ ,$$for some polynomial $J_{m,i}(x)$. By first applyingLemma~\ref{lemmf.1}for $u(x)=D_{m,i+1}(x)$, $v(x)=J_{m,i}(x)$,$s(x)=D^2_{m,i}(x)G_{m+1}^*(x) $ and thenapplying Fact~\ref{factA.2} for $q(x)=F_{m+1}^*(x)$,$r(x)=G_{m+1}^*(x)$, we obtain that\begin{eqnarray}\label{equaf.29}\|D_{m,i+1}(x)\|\leq(256/255)\|D_{m,i}(x)\|^2(k/\eta)   \|G_{m+1}^*(x)\|\ \|F_{m+1}^*(x)\|\nonumber\\\leq (257/255)2^{n-1}\|D_{m,i}(x)\|^2k/\eta.\end{eqnarray}Due to a sufficiently small initial upper bounds of (\ref{equaf.24}) on the norm of $D_{m,0}(x)$  and of (\ref{equaf.5}) on $\delta_0$ the bound (\ref{equaf.29}) implies a rapid decrease of the norm $\|D_{m,i}(x)\|$with the growth of $i$, so that already for $i\leq 3$, we will arrive at (\ref{equaf.3})--(\ref{equaf.5}) with $m$ replaced by $m+1$, $H_{m+1}(x)=H_{m,i}(x)$, and $D_{m+1}(x)=D_{m,i}(x)$.Indeed, our only remaining task is to ensure that$\|D_{m,3}(x)\|\leq \delta_m^{1.5}$, and due to(\ref{equaf.29}), this only requires to satisfy the bounds$$ \|D_{m,3}(x)\|\leq (257/255)(k/\eta)2^{n-1}\|D_{m,2}(x)\|^2\leq (257/255)^3(k/\eta)^32^{3n-3}\|D_{m,1}(x)\|^4\leq $$$$(257/255)^7(k/\eta)^72^{7n-7}\|D_{m,0}(x)\|^8\leq\delta_m^{1.5}.$$Due to (\ref{equaf.24}), the latter requirement issatisfied if $$(257/255)^7(k/\eta)^72^{7n-7}(2.1)^8\delta_m^82^{8n+8k}(k/\eta)^{16}\leq \delta_m^{1.5} $$or, equivalently, if $$1/\delta_m^{13/2}\geq (257/255)^7(2.1)^8(k/\eta)^{23}2^{15n+8k-7}.$$Since $(257/255)^{14/13}(2.1)^{16/13}<2.52$, it sufficesto choose any $\delta_0$ satisfying the bound $$1/\delta_0\geq 2.52(k/\eta)^{46/13}2^{(30n+16k-14)/13},$$which is even substantially milder than the bound of(\ref{equaf.5}) on $\delta_0$. This completes the description of part b) of thecomputation, whose arithmetic cost is within thebound of Proposition~\ref{propA.1}, as simple inspection shows.Moreover, the introduction of the additional iterativeloop of computing $H_{m,i}(x)$, $D_{m,i}(x)$ enabled us toweaken the initially assumed bound (\ref{equaf.5}) on $\delta_0$. Even milder bounds on $\delta_0$would suffice if we choose to stop in $i$ steps of the iteration(\ref{equaf.28a}), (\ref{equaf.28}) for $i>3$, setting $H_{m+1}(x)=H_{m,i}(x)$, $D_{m+1}(x)=D_{m,i}(x)$. Furthermore, simpleinspection shows that,unlike the algorithm of Appendix A, the approachof this appendix does not require to increase theprecision of computation above $O(bn)$ bits.Moreover, at some stages of computations even a lowerprecision suffices. \newpage \section{Figures.}\label{secg}\setcounter{equation}{0} \begin{figure}[h]\begin{center}\setlength{\unitlength}{1in}\begin{picture}(5,5)\multiput(0,0)(5,0){2}{\line(0,1){5}}\multiput(0,0)(0,5){2}{\line(1,0){5}}\multiput(1.5,1.5)(2,0){2}{\line(0,1){2}}\multiput(1.5,1.5)(0,2){2}{\line(1,0){2}}\put(2.5,2.5){\vector(1,0){1}}\put(2.5,2.5){\vector(0,1){2.5}}\put(2.57,4.2){R}\put(3,2.4){r}\put(2.45,2.37){X}\multiput(1.7,1.74)(0.23,0.25){5}{\circle*{0.015}}\multiput(3.3,3.25)(-0.05,-0.3){5}{\circle*{0.015}}\end{picture}\caption{\label{figu1} Isolation  ratio for squares,$i.r.(S(X,r)).$}The dots show all the zeros of $p(x)$ lying in the larger square.They also lie in the smaller square. \end{center}\end{figure} \newpage \begin{figure}[h]\begin{center}\setlength{\unitlength}{1.25in}\begin{picture}(6,6)\put(0,0){\framebox(4,4){}}\put(0,2){\line(1,0){4}}\put(2,0){\line(0,1){4}}\multiput(1,2)(2,-2){2}{\line(0,1){2}}\multiput(2,1)(-2,2){2}{\line(1,0){2}}\multiput(1,3.5)(1,-3){2}{\line(1,0){1}}\multiput(1.5,3.0)(1,-3){2}{\line(0,1){1}}\multiput(1.17,3.15)(0.23,0.25){2}{\circle*{0.01}}\multiput(2.45,0.35)(0.25,0.3){2}{\circle*{0.01}}\end{picture}\caption{\label{figu2}Weyl's algorithm.}The dots show all the zeros of $p(x)$ lying in the large square. \end{center}\end{figure}  \newpage \begin{figure}[h]\begin{center}\setlength{\unitlength}{0.8in}\begin{picture}(6.4,6.4)%\put(0,0){\framebox(6.4,6.4){}}\multiput(0,0)(6.4,3.2){2}{\line(0,1){3.2}}\multiput(0,0)(3.2,6.4){2}{\line(1,0){3.2}}\multiput(1.6,0)(3.2,1.6){2}{\line(0,1){4.8}}\multiput(0,1.6)(1.6,3.2){2}{\line(1,0){4.8}}\put(0,3.2){\line(1,0){6.4}}\put(3.2,0){\line(0,1){6.4}}\multiput(1.6,1.2)(.4,.4){10}{\line(0,1){0.4}}\multiput(1.2,1.2)(.4,.4){10}{\line(1,0){0.4}}\multiput(4.8,5.2)(-0.4,-0.4){10}{\line(0,-1){0.4}}\multiput(5.2,5.2)(-0.4,-0.4){10}{\line(-1,0){0.4}}%\multiput(2,1.2)(1.6,1.6){3}{\li\multiput(2,1.2)(1.6,1.6){3}{\line(0,1){0.4}}\multiput(1.6,1.2)(1.6,1.6){3}{\line(1,0){0.4}}\multiput(1.2,2)(1.6,1.6){3}{\line(0,-1){0.4}}\multiput(1.6,2)(1.6,1.6){3}{\line(-1,0){0.4}}\end{picture}\caption{\label{figu3}16 smaller suspect squares at}Stage $g+2$ versus 10 larger ones \\at Stage $g$ of Weyl's algorithm.\end{center}\end{figure}  %%\newpage%%\begin{figure}[h]%%\begin{center}%%\setlength{\unitlength}{1.25in}%%\begin{picture}(6,6)%%%%%%\put(3,1){\line(2,-1){2}}%%\put(3.78,0.5){$\rho$}%%%%\put(0,1){\line(0,1){3}}%%\put(0,1){\line(1,0){5}}%%\put(5,0){\line(0,1){1}}%%\put(0,4){\line(3,-3){3}}%%\put(0,4){\line(5,-3){5}}%%%\put(5,0){\line(-2,1){2}}%%\put(5.3,-0.5){$\geq \frac{L}{2}$}%%\put(0,4.08){$\beta_h$}%%\put(3,0.9){$\alpha_h$}%%\put(5,-0.04){$c$}%%\put(4,1.05){$\geq \frac{L}{2}$}%%\put(-0.6,2.5){$2z_h/\sqrt 2$}%%\put(1.5,0.8){$2z_h/\sqrt 2$}%%%\put(4,0.4){$r$}%%\put(1.5,2.5){$2z_h$}%%\put(3.5,1.9){$\Delta_h$}%%%%\end{picture}%%\caption{\label{figu5} Proof of Claim~E.1: Projecting}%%the segments $[\alpha_h, \beta_h]$ and $[\alpha_h, c]$%%\\%%onto the real axis.%%\end{center}%%\end{figure} \end{document}