From:	IN%"sad@cunyvms1.gc.cuny.edu" 29-FEB-1996 14:08:56.64To:	IN%"vpan@lcvax.lehman.cuny.edu"CC:	Subj:	This is the latest file on Parallel Comp. of poly. GCD and rel. \\ comp.over abstract field.Return-path: <sad@cunyvms1.gc.cuny.edu>Received: from vice1.lehman.cuny.edu by lcvax.lehman.cuny.edu (PMDF #12809) id <01I1SDY2WFRK0015KW@lcvax.lehman.cuny.edu>; Thu, 29 Feb 1996 14:08 ESTReceived: from CUNYVMS1.GC.CUNY.EDU by vice1.lehman.cuny.edu (5.51/SMI-3.0DEV3) id AA03376; Thu, 29 Feb 96 13:54:08 ESTReceived: by cunyvms1.gc.cuny.edu (MX V3.3 VAX) id 4401; Thu, 29 Feb 1996 13:45:46 ESTDate: Thu, 29 Feb 1996 13:44:56 ESTFrom: sad@cunyvms1.gc.cuny.eduSubject: This is the latest file on Parallel Comp. of poly. GCD and rel. \\ comp.over abstract field.To: vpan@lcvax.lehman.cuny.eduMessage-id: <0099EA34.936E4BCC.4401@cunyvms1.gc.cuny.edu>\documentstyle [11pt] {article}\topmargin-.25in\oddsidemargin.05in\evensidemargin.05in\textheight8.55in\textwidth6in\renewcommand{\baselinestretch}{1.5}\begin{document}\begin{center}{\bf Parallel Computation of Polynomial GCD and Some Related Parallel Computations over Abstract Fields}\end{center} \begin{center}{\bf  Victor~Y.~Pan \\Department of Mathematics and Computer Science\\ Lehman College, City University of New York\\ Bronx, NY 10468 \\Internet:  VPAN@LCVAX.LEHMAN.CUNY.EDU}\end{center}\thispagestyle{empty}\begin{center}{\bf Summary  \footnote{The results of this paper have been presented at the 7th Annual ACM-SIAM Symposium on Discrete Algorithms, in Atlanta, Georgia, in January 1996; the author is greatful to ACM and SIAM for permission to cite them.}}\end{center}%(over any field of constants) Several fundamental problems of computations with polynomials and structured  matrices are well known for their resistance to effective parallel solution. This includes computing the gcd, the lcm, and the resultant of two polynomials, as well as any selected entry of either the extended Euclidean scheme for these polynomials or the Pad\'{e} approximation table associated to a fixed power series, the solution of theBerlekamp-Massey problem of recovering the $n$ coefficientsof a linear recurrence sequence from its first $2n$ terms,the solution of a Toeplitz linear system of equations, computing the ranks and the null-spaces  of Toeplitz matrices, and several other  computations with Toeplitz, Hankel, Vandermonde, and Cauchy (generalized Hilbert) matrices and with matrices havingsimilar structures. Our algorithms enable us to reach new record estimates for randomized parallel arithmetic complexity of all these computations (over any field of constants), that is, $O((\log n)^h )$ timeand $O((n^2 \log\log n)/(\log n)^{h-g} )$ arithmetic processors, where $n$ is the input size,  $g$ varies from 1 to 2,  and $h$ varies from 2 to 4, depending on the computational problem and on the fixed field of constants. The estimates include the cost of correctness tests for the solution. The results have further applications to many related computations over abstract fields and  rely on some new techniques (in particular, on recursive Pad\'{e} approximation and regularization of Pad\'{e} approximation via randomization), which might be  of some independent technical interest. In particular, an auxiliary algorithm computes (over any field) the coefficients of a polynomial from the power sums of its zeros, which is an important problem of algebraic coding.%\end{abstract} \noindent{\bf Key Words:}Toeplitz linear systems, polynomial gcd,Toeplitz-like matrices, Newton's identities,finite fields, randomized algorithms, algebraic computations, parallelcomputational complexity.\vspace{.23in} \noindent{\bf  1991 Mathematics Subject Classification:}68Q40, 68Q25, 15A06, 68Q22%\newpage\vspace{.3in} \noindent{\bf 1. Introduction.}\vspace{.2in} \noindent {\bf 1.1 Our subject.} Fundamental computations with Toeplitz, Hankel, Toeplitz-like, and Hankel-like matrices (generalizing Toeplitz and Hankel matrices) as well as with some other dense and structured matrices are among the major topics of computations in linear algebra (for definitions, see section 2 or [P89/90] and [BP94]). These computations include solving nonsingular and singular linear systems of equations with dense and structured coefficient matrices $V$,  computing the ranks of such matrices, their null-spaces, determinants, and inverses (where $\det V \neq 0$). The algorithms for these  problems have numerous extensions and applications. In particular, in section 12, we consider extensions to the celebrated problems of computing the greatest common divisor (gcd), the least common multiple (lcm), the resultant, and all the entries of the extended Euclidean scheme for a pair of univariate polynomials, of computing a fixed entry  of the Pad\'{e} table for a given formal power series, and of recovering the $n$ coefficients  of a linear recurrence sequence from its first $2n$ terms (known as the Berlekamp-Massey problem). In sections 2 and 12, we will recall some reductions of these problems to Toeplitz and/or Toeplitz-like computations and their     correlations to each other; this material can be found in [BGY80], [BGH89], [G84], [BP94] as well (also compare [Col67], [Bro71]). In particular, facts 12.1 and 12.3 of section 12 enable us to reduce the gcd problem at first to solving a simple polynomial identity and then, further, to Toeplitz-like computations, namely, to computing the rank of a Sylvester (resultant) matrix and to solving a nonsingular (resultant or subresultant) linear system of equations. Alternatively, the gcd problem can be easily reduced to computing an entry of a Pad\'{e} table for the formal power series defined by the quotients of 2 polynomials $u(x)/v(x)$ or $u(1/x)/v(1/x)$ (see our appendix A), and this problem is equivalent to the rank computation and linear system solving problems for Toeplitz input matrices (of [BGY80] or [BP94]). Due to these and other known reductions, we will focus on Toeplitz and Toeplitz-like computations, in particular, on the rank computation and linear  system solving. We will conclude this subsection with recalling about some further major applications of Toeplitz-Hankel type computations to the theory  and practice of computing for sciences, engineering, and communication. They include applications to  shift register synthesis, BCH and Read-Solomon decoding of error-correcting codes, inverse scattering, adaptive filtering, modelling of stationary and nonstationary processes, numerical computations for Markov chains, sparse multivariate polynomial interpolation, solution of partial differential and integral equations, parallel computations with general matrices over an arbitrary field of constants, and approximating polynomial zeros [Be68], [Bro71], [BT71], [Col67], [K74], [Ma75],  [KVM78], [BGY80], [L-AK84], [L-AKC84], [BK87], [K87], [KP91], [P91], [P92a], [KP92], [PSLT93], [BP94], [Ste94], [P95], [P,b]. \vspace{.3in} {\bf 1.2. Model of Computation.} We will be looking for fast and processor efficient parallel algorithms for the listed computational problems. Stating our complexity estimates, we will assume the customary EREW PRAM model of parallel computing [EG88], [KR90], [J92], under which, at every time-step, each non-idle processor may perform a single unit cost operation (in our case, arithmetic operation or a comparison with 0). In fact, to obtain our upper estimates for the cost of performing our algorithms over the complex field of constants, we essentially need just a model that supports performing fast Fourier transform at the $n$-th  roots of 1 by using $O(\log n)$ time on $n$ processors, and this can be achieved, for instance, on the butterfly, shuffle-exchange, or hypercube processor arrays [Le92], [Q94].  Over an arbitrary field of constants, our cornerstone operation will be polynomial multiplication modulo $x^n$, which can be performed by using $O(\log n)$ time and $O(n \log \log n)$ processors (cf. [CK91]) under various realistic computer models too. Hereafter, we will write $O(t,p)$ in order to denote simultaneous upper bounds $O(t)$ on time and $O(p)$ on the number of processors used. We will assume a variant of {\it Brent's scheduling principle} [Bre74], according to which a single processor may simulate one step of $s$ processors in $O(s)$ steps, so that $O(t,p)$ implies $O(tp,1)$, but not necessarily vice versa. The product $w=tp$ of parallel time and the number of processors used by a parallel algorithm is called its {\it potential work} [KR90], [PP95]. The algorithm is called {\it work (or processor) efficient} if it supports the bound $O((\log n)^d)$ on the {\it work inefficiency ratio} $w(n)/t_{seq}(n) $, where $n$ is the input size, $w(n)$ is the work of this algorithm, $t_{seq}(n)$ is the current record sequential time for the solution of the same problem, and $d$ is a constant, independent of $n$. If the ratio remains $O((\log n)^d)$ even when we replace the current record $t_{seq}(n)$ by a known lower bound on the sequential time, the algorithm is called {\it strongly work (or processor) efficient}.  An algorithm is {\it in NC} if its parallel time is $O((\log n)^d)$ and if its work and processor bounds are $O(n^c)$ for 2 constants $c$ and $d$, independent of $n$.  {\it RNC} stands for {\it randomized NC}; the randomized estimates for the computational cost are of {\it Las Vegas type}, if they include the cost of testing correctness of the output, and of {\it Monte Carlo type}, otherwise. Devising  algorithms that would be both in NC or RNC and work efficient turned out to be a quite difficult task for the computational problems cited in the previous subsection.\vspace{.3in} \noindent {\bf 1.3. Fast sequential solution and its parallel acceleration based on algorithms for computations with general matrices.} The sequential time bound $$t_{seq}(n) = O((\log n)^2 n \log \log n ) \eqno (1.1)$$is obtained (over any field of constants) for the gcd, lcm, Pad\'{e} and Berlekamp-Massey problems of size $n$, due to their solution by means of a fast version of the Euclidean algorithm [AHU74], [BGY80], where the extra factor $\log \log n$ is the (small) price for using fast polynomial multiplication over any field of constants [CK91]. The same sequential time bound, for solving a nonsingular Toeplitz or Toeplitz-like linear system of $n$ equations, has been supported by the algorithms of [BGY80] for a Toeplitz system and of [K94], [K,a] for a Toeplitz-like system, the latter algorithm relies on [BA80], [Morf80] (compare also [Mus81], [dH87], [AGr88], [ChK91], [GG94]). All these fast sequential algorithms, however, {\it resist their effective parallelization.} Moreover, for some time, the only known way to obtaining NC algorithms for the above computational problems was by first reducing each of these problems to solving $O(n)$ nonsingular Toeplitz or Toeplitz-like systems, each of $n$ equations [hereafter referred to as $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$] (compare sections 2 and 12 and [BGH82], [G84]) and then by applying the known fast parallel algorithms available for solving a general nonsingular linear system of equations [Cs76], [PSa78], [Be84], [Ch85], [GP89]. The latter algorithms support deterministic cost bounds $O((\log n)^2, n^{\omega + \sigma})$, where $\omega \geq 2$ is the exponent of matrix multiplication, whose current record upper bound is $\omega < 2.375...$ and where $\sigma$ varies from 1/2 [PSa78] or slightly less than 1/2 [GP89] (over fields of characteristics $p=0$ and $p>n$) to 1 (over any field [Be84], [Ch85]). The time bound $O((\log n)^2)$ is quite attractive, but the processor and work bounds of order $n^{\omega + \sigma}$ for a single linear system are clearly too high, in comparison with (1.1).  \vspace{.3in} {\bf 1.4. A solution over the fields of characteristics 0 and greater than $n$.} Over the fields of characteristics $p=0$ and $p>n$, the algorithm of [P90/92] yields the bounds $$ C_n = O((\log n)^2, (n^2 \log \log n)/\log n)  \eqno (1.2)$$ on the parallel cost and $ w(n)/t_{seq}(n) = O(n)$ on the work efficiency ratio  for  $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$. For an $n \times n$ Toeplitz-like matrix $V$ and a vector $\vec{f}$, the algorithm computes the Krylov sequence $\{V^i \vec{f}\}$ and the values $\{ {\rm trace} \ (V^i) \}$, $i=0, 1, \ldots, m$, $m=O(n)$, at the cost $C_n$ of (1.2), over any field of constants. Given the traces of $V^i$, a simple and elegant algorithm of [Sc82] (see our section 3) computes the coefficients of $C_V(x)$, the characteristic polynomial of $V$, at the cost $O((\log n)^2, (n^2 \log \log n)/\log n)$, over the fields of characteristics $p=0$ and $p>n$. Given $C_V(x)$,  one easily obtains rank $V$ (if $p=0$), $\det V$, and if $\det V \neq 0$, then also $V^{-1}$ and $V^{-1}\vec{f}$, all at the cost $C_n$ of (1.2). If $p=0$, then, within the same deterministic cost bound $C_n$ of (1.2), the algorithm of [P90/92] can be extended to compute the rank of $V$, via computing the characteristic polynomial of the Toeplitz-like symmetric matrix $\left( \begin{array}{lr} O & V \\ V^T & O  \end{array} \right)$. If $p>n$, we may extend the algorithm of [P90/92], so as to compute the rank at a Las Vegas randomized cost $ \hat{C}_n=O((\log n)^2, (n^2 \log \log n)) %\ .  \eqno (1.3)$.  Then, we may easily solve the gcd, lcm, Pad\'{e}, and Berlekamp-Massey problems at the latter (randomized) cost for $p>n$ and at a deterministic cost bounded by (1.2) for $p=0$. The  processor bound % and by (1.3) for $p>n$. The processor bound of $\hat{C}_n$ decreases by factor $\log n$, over the fields containing sufficiently many distinct elements, and similarly in our randomized cost bounds (1.3), (1.4), and (1.5) of subsections 1.5 and 1.6.\vspace{.3in} {\bf 1.5. Extensions to computations over an arbitrary field (previous works).}  The results of [P90/92] do not apply, however, in the case of positive characteristics $p \leq n$, due to the problem of the transition from $\{ {\rm trace} \ (V^i) \}$ to $C_V(x)$. (Note that $p \neq 1$, so $p\geq 2$ if it is positive.) This problem [which we call the {\it inverse power sum }  ($ \hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n)$) problem] has also been well known as a bottleneck of some important algebraic algorithms for coding ([Be68], pp. 178--179, cites  $ \hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n)$ for $p=2$) as well as of general matrix computations, where this problem restricts, for instance, the algorithms of [Cs76], [PSa78], [GP89], and [KP91] to the cases of $p=0$ and $p>n$.  Two approaches proposed for the solution of the  $ \hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n)$ problem over any field of constants were based on the techniques of {\it recursive triangulation } [KP92] and of finding {\it fundamental sets of power sums} [V900], [K25], [K27], [N27], [Sc93], [P93a], respectively. (Besides, both approaches also relied on some techniques of {\it matrix regularization}.) For general linear systems, both approaches yielded NC and processor efficient randomized algorithms over any field (see [KP92] on the former approach and [P93a] or [BP94], appendix C of chapter 4, on the latter approach). In [P93a] and [BP94] (proposition 4.6.4 of chapter 4), the technique of recursive triangulation has also been proposed in order to prove the Las Vegas randomized bound $$\tilde{C}_n = O((\log n)^3, (n^2 \log \log n)/\log n) \eqno (1.3) $$ on the complexity of $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$ over any field of constants. The proof, however, has an error; a nontrivial modification of this approach in [KP94]  (still based on recursive triangulation) has salvaged the bound (1.3) for computations over the fields whose characteristics $p$ satisfy the inequality $p \leq n^{\frac{1}{\omega -1}}$, $\omega \geq 2$ being the exponent of matrix multiplication. Since the best available upper bound on $\omega$ is $2.375\ldots$ [CW90], the latter result of [KP94] has still left the problem of proving (1.3) open in the cases where $n^{0.72\ldots} =n^{\frac{1}{2.375\ldots}}<p \leq n$. \vspace{.3in} {\bf 1.6. Our results and organization of the paper.} In the present paper, we solve the latter problem based on the extension to the Toeplitz-like case of the version of the fundamental set techniques proposed for general matrix computations in [Sc93], [P93a], and appendix C of chapter 4 of [BP94]. We organize our presentation as follows. After some preliminaries in section 2, we recall the solution of the $ \hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n)$ problem, for $p=0$ and $p>n$, in section 3. In sections 4 and 5, we extend this solution to the case of any $p$. (Recall the cited applications to algebraic codes in [Be68] in the case $p=2$ and see our remarks 5.3 and 7.3 in sections 5 and 7.) Consequently (due to the algorithm of [P90/92] and to some new techniques of {\it recursive Pad\'{e} approximation} and of {\it regularization of Pad\'{e} approximation via randomization}), we arrive at the Las Vegas randomized complexity bounds (1.3) for $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$ over any field. In section 6, we easily extend these complexity bounds to the problem of computing  $\det V$, the determinant of a Toeplitz-like matrix $V$, over any field, and if $\det V \neq 0$, then also to the problem of computing $V^{-1}$. In section 7, we elaborate the extension of (1.3) to the computation of rank $V$. Although such an extension has been claimed in [P93a] and [BP94], the proofs have some gaps, which we fill, by means of extending some special randomization techniques  of [P93a] (also presented in [BP94])  and applying {\it double loop} of our techniques for  recursive Pad\'{e} approximation (with its regularization via randomization). Due to doubling the loop, our parallel time bound for Toeplitz and Toeplitz-like rank computation increases by factor $(\log n)/\log p$ for $2\leq p \leq n$, versus the solution of $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$, so that over any field, we arrive at the Las Vegas randomized cost bound $$\hat{C}_n=O((\log n)^4, (n^2\log \log n)/(\log n)^2) \ . \eqno (1.4)$$  In sections 8 and 9, we extend the results of sections 5, 6, and 7 in order to compute a basis for the null-space, ${\bf N}$(V), of a Toeplitz-like matrix $V$ and a specific solution $\vec{x}$ to a singular Toeplitz-like linear system of equations $V\vec{v}=\vec{f}$, whose general solution set can be represented as $\vec{x}+ {\bf N}$(V). The cost bounds, (1.2) for $p=0$, (1.2) or (1.3) for $p>n$, and (1.4) for any $p$, are preserved, though randomization is required even for $p=0$ in this case. In section 10, we extend these results to computation of the Moore-Penrose generalized (pseudo) inverse $V^+$ of a Toeplitz-like matrix $V$ (we show that $V^+$ is a Toeplitz-like matrix too) and of a least-squares solution $\vec{x}=V^+\vec{f}$ to a singular Toeplitz-like linear system $V\vec{x}=\vec{f}$. (Defining $V^+$ and the least-squares solutions, we restrict our study in section 10 to the computations in the field of complex numbers ${\bf C}$ or in its subfield.)  In section 11, we recall that the techniques of [P89/90] and [BP93] enable us to extend our results to similar fundamental computations for some other major classes of structured matrices (Cauchy-like, Vandermonde-like, and Toeplitz-like + Hankel-like). The elaboration of the impact on the computations with general matrices is postponed, to be presented in [KP,a] and/or [BP,a]. In section 12, the deterministic bound (1.2) for $p=0$ and the randomized bounds (1.3) and (1.4), for any $p$, are extended from $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$ and Toeplitz-like rank computation to the gcd, lcm, resultant, Pad\'{e} and Berlekamp-Massey problems, and we also compute the extended Euclidean scheme for 2 polynomials of degrees at most $n$ at deterministic cost $O((\log n)^3, (n^2\log \log n)/(\log n)^2)$ over the fields of characteristic $p=0$, at randomized cost (1.4) for $p>n$, and at randomized cost $$C_n^* = O((\log n)^5, (n^2\log \log n)/(\log n)^3) \  \  {\rm for} \  0< p \leq n \ . \eqno (1.5)$$ \vspace{.3in} {\bf 1.7. Discussion.}The presented NC and RNC algorithms enabled us to decrease processor inefficiency of the known NC and RNC algorithms for some major computational problems over abstract fields. The extended Euclidean scheme and also the inverse of a Toeplitz or Toeplitz-like matrix, may have order $n^2$ distinct entries. Consequently, the sequential time bounds for their computation are at least of order $n^2$, which implies strong work (or processor) efficiency of our parallel NC and RNC algorithms in these cases. For other computational problems, we yielded work bounds of order $n^2$, versus a known sequential time bound of order $n$ (in both cases, up to polylogarithmic factors); devising NC  or RNC and work efficient algorithms for these problems (over abstract fields) remains a research challenge. An encouraging recent result is the design of RNC and work efficient algorithms for these problems in the field of rational numbers provided that rouding-off a rational number to the closest integer is allowed as a unit cost operation. In this case, a nontrivial parallel algorithm of [P,a], computes the gcd of two polynomials with rational coefficients and solves the other cited computational problems at the randomized parallel arithmetic cost  $$c_A(n)=O((\log n) \log a, n \log n)  \ , \eqno (1.6)$$by using arithmetic operations with rational numbers, each represented with  $O(n\log a)$ bits, provided that the input is represented with $O(n)$ parameters and with $\log a$ bits, $a \geq n$. The algorithm of [P,a] relies on a combination of the variable diagonal techniques of [P85], [P87] as well as some other techniques of [P85], [P87] (proposed for computations with general matrices) with some more recent advanced techniques of [P92c], [P93], and [P93b] for computations with Toeplitz and Toeplitz-like matrices.  Technically, [P,a] extends and elaborates a similar algorithm of [BP94], page 357, which supports randomized cost bounds $O((\log n)^2\log a, n)$, [slightly inferior to (1.6)] and the same bit-precision bounds. In [R95], the cost bounds  $O((\log n)^2, (\log n)^{2.376} n)$ have been claimed, under the restriction that all entries of the input Toeplitz-like matrices are integers in the range from $-n^d$ to $n^d$ for a fixed constant $d$, independent of $n$. The presentation, statements of the results, and historical account are quite inaccurate in [R95]. For instance, randomization is used in the main algorithm of [R95] but not referred to in the statements of the complexity estimates. In another instance, the upper estimates of lemma 2.4 on page 753 are fundamental for the main algorithm of [R95], but their proofs are absent from both [R95] and the bibliography cited by Reif in conjunction with lemma 2.4 (in fact, part of this bibliography is irrelevant). Furthermore, the upper estimates stated in the part b) of the lemma are much worse (much greater) than the known ones, whereas the upper estimates stated in the part c) are '' better'' (less) than the best record estimates presently available.  There are even more major problems with the paper [R95], however, which strongly suggest that the latter  paper had actually presented the results not yet obtained by the author by the moment of its appearance: \\ a) the algorithm description and proofs are a lengthy repetition of the previous works and known techniques (particularly, of ones from [P85], [P87], [P88], [P90], [P90/92], and [PR91], though the sources are not properly cited) but with wide gaps and omissions whenever new results are needed; \\b) there are most serious technical errors in the attempted substantiation of the main claimed results.In appendix B, we comment on one group of such major errors and flaws;  this group actually implies the increase by factor $n^{\omega -1}$ of Reif's claimed bounds on work-complexity, where $\omega$ is the exponent of matrix multiplication, $2\leq \omega < 2.376$. Another major area of technical failure of [R95] is control over the bit-precision of the computations. Most of the necessary estimates are simply absent from [R95], and what is presented, again, has most serious technical errors. In particular, due to one such technical error (see our specific comments in [P,a]), the bit-precision of the computations in [R95] and, consequently, their Boolean complexity bounds are actually by factor $n$ greater than they are claimed to be. Moreover, fixing the latter error of [R95] (as well as fixing the errors of [R95] pointed out in appendix B of our present paper) is nontrivial and requires to use techniques that are substantially more sophisticated than any technique proposed in [R95]. Incidently, [R95] cites [R93], which has similar problems: highly inaccurate presentation, most serious errors in the proofs, and technically wrong statements; moreover, the main claimed result (which is actually the only result of [R93]) is a rediscovery of a result from the earlier publications [B-OT90] and [BP91] (where, unlike [R93], the proofs are correct). Besides arithmetic operations, the algorithms of [P,a], [BP94], and [R95] use rounding to the closest integer, and this does not allow one to apply the powerful techniques available for algebraic and symbolic computations  [such as (Newton-Hensel's)  $p$-adic lifting and/or the Chinese remainder algorithm] and customarily applied  in order to decrease the precision of computing. To be able to use such techniques,  we would have had to restrict the model of computing by allowing only arithmetic  (field) operations, but then the known record bounds on the parallel arithmetic complexity of the solution increase to the level of (1.2)--(1.5).  Of course, besides the goal of reaching an NC and processor efficient solution, even more minor progress would be interesting. For instance, the present paper raises the question if double recursion for Pad\'{e} approximation can be avoided; using single recursion would have enabled us to decrease the time bound by factor $(\log n)/\log p$. Another open problem suggested by this paper is to prove some of the bounds (1.3)--(1.5) as deterministic complexity bounds  or at least to decrease the number of random parameters involved in the algorithms supporting these bounds and the failure probability of these randomized algorithms.   \vspace{.25in} \noindent{\bf 2. Auxiliary results and definitions.} All our randomized complexityestimates do not include the cost of generating random parameters,base 2 will be assumed for all logarithms, and $\log^*n$ is definedas follows: $\log^{(0)} n =n, \log^{(h+1)}n=\log \log^{(h)}n,\log^*n=max\{h:\log^{(h)}n>0\}$. We will use some known results and customarydefinitions, which can be found in[BP94]. {\bf Fact 2.1} [CK91]. {\it The coefficients of the product and the sum of two univariate polynomials of degrees at most $n$ can be computed over any field of constants at the cost bounded by $O (\log n, n \log\log n)$, for the product,and $O(1,n)$, for the sum.} {\bf Fact 2.2} [BP 92/93]. {\it The coefficients of a polynomial $(1/B(x)) \  {\rm mod} \ x^n$, where $B(x)$ is a given univariate polynomial of a degree $n$ such that$B(0) \neq 0$, can be computed over any field of constants at the cost bounded by$O((\log n) \log^* n, (n \log\log n)/\log^* n)$ or,alternatively, by $O(\log n, (\log\log n)^2 n)$.} {\bf Remark 2.1} [CK91]. The estimates of facts 2.1 and 2.2 hold  for $k$-variate polynomials of degrees at most  $n_i$ in the $i$-th variable, $i=1, \ldots, k$, for $n=\prod_{i=1}^k n_i$, and for $x$ replaced by $\vec{x} =(x_1, \ldots, x_k)$.  {\bf Definition 2.1.} An $n\times n$ matrix $ T=(t_{i,j})$ is a{\it Toeplitzmatrix} if $t_{i+1,j+1} = t_{i,j}$ for all pairs $i,j$, where$0 \leq i,j \leq n-1$.  $ L ( \vec{z} ) $ is the $ n \times n $ lowertriangular Toeplitz matrix with the first column $ \vec{z}$. $W^T$ is the transpose of a matrix or a vector $W$. $Z$ is the triangular Toeplitz matrix $L(\vec{e}^{\ (1)})$,  where $\vec{e}^{\ (1)} = (0, 1, 0, \ldots, 0)^T$;$O$ is a null matrix of an appropriate size; $I= \displaystyle{ \left( \begin{array}{ccc}1 & & O \\ & \ddots & \\ O & & 1 \end{array} \right) }$and $ J = \displaystyle{ \left ( \begin{array}{lcccr}  O  &  &  &  & 1 \\ &         &         & {\bf .} &    \\ &         & {\bf .} &         &  \\  &  {\bf .}&  &  &    \\  1     &  &  &  &  O                   \end{array} \right ) }$ are the$n\times n$ identity and reversion matrices, respectively, so that $Z\vec{v} = (0, v_1, \ldots, v_{n-1})^T$, $I \vec{v} = \vec{v}$,$J \vec{v} = (v_n ,\ldots ,v_1 )^T$ for any vector$\vec{v} = (v_1 ,\ldots ,v_n )^T$, $J^2=I$.$H$ is called a {\it Hankel matrix} if $JH$ is a Toeplitz matrix. {\bf Definition 2.2.} An $n \times n$ matrix $V$ is {\it Toeplitz-like} if itis given as a sum $ \sum_{i=1}^{\ell} L ( \vec{x} _i )L^T ( \vec{y}_i ) $ of $\ell = O(1)$ products $ L ( \vec{x}_i )L ^T ( \vec{y}_i ) $ of pairs of triangular Toeplitzmatrices or, equivalently (cf. [KKM79]), if $V-ZVZ^T = \sum_{i=1}^{\ell} \vec{x}_i \vec{y}_i^T$.The pair of $n \times \ell$ matrices $ X = ( \vec{x}_1 , \ldots , \vec{x}_{\ell} ) ,Y = ( \vec{y}_1 , \ldots , \vec{y}_{\ell} ) $ is called a {\it displacement generator for $V$ of length } $ \ell .$ The minimum length $\rho=\rho_V$ of a displacement generator for a matrix $V$ is called the {\it displacement rank} of $V$. (It is possible to decrease the length $\ell$ of a given displacement generator of an $n \times n$ matrix $V$ to the level $\rho$, at a cost $O(\log n, n\log \log n)$ if $\ell=O(1)$, see appendix B.)A matrix $JV$ is called {\it Hankel-like} if $V$ is a Toeplitz-like matrix. Clearly, any Toeplitz matrix can be immediately represented with  its displacement generator of a length at most 2, and short displacement generators can be obtained for product and sums of a few Toeplitz matrices, as well as for their inverses (see fact 2.4 below). Many properties of Toeplitz matrices can be extended to Toeplitz-like matrices. In particular, representing Toeplitz-like matrices with their displacement generators, rather than with their entries, saves the memory space, whereas operating with such generators,  instead of operating with matrices, saves computational time and/or processors to be involved. In this paper, we will need only a small part of the wealth of the techniques available for computations with Toeplitz-like and other structured matrices. Many results and techniques known in this area can be found in [BP94], and we also refer the reader to [KKM79], [BA80], [Morf80], [CKL-A87], [KP91], [KP92], [P91], [P90/92], [P92a], [P93], [P93b], and [BP93] on such techniques, further results, and applications.{\bf Remark 2.2.} There are also other classes of displacement generators (with similar properties) for Toeplitz-like matrices, for Hankel-like matrices, and for Toeplitz-like + Hankel-like matrices, each represented as a sum of 2 matrices: of a Toeplitz-like one and of a Hankel-like one (see [KKM79], [BP93], [BP94]). Likewise, appropriate scaling operators and displacement and scaling operators are naturally associated with Cauchy (generalized Hilbert) and Vandermonde matrices, respectively, and such an association is immediately extended to the larger classes of Vandermonde-like and Cauchy-like matrices represented by their short scaling and displacement-scaling generators, respectively (compare [P89/90], [BP94]). Such representations enable economical storage and effective computations for the structured matrices of all these classes. {\bf Fact 2.3.} {\it Any $n\times n$ Toeplitz-like matrix $V$ can be multiplied by a vector $\vec{u}$ at the computational cost$O(\log n,n \log\log n)$, over any field of constants.}{\bf Proof.} Observe that $V\vec{u}$ is a subvector of the coefficient vector of the product of 2 polynomials one of whose coefficient vectors is just the vector $\vec{u}$ and another is  formed by the entries of the matrix $V$ ([BP94], page 133). Then, apply fact 2.1. \hfill  $\Box$  \vskip2mm{\bf Fact 2.4.} {\it Let $T$ and $V$ be a pair of $ n \times n $Toeplitz-like matrices given with their displacement generators of lengths $ \ell _T $ and $ \ell _V ,$  respectively.Then the matrix $TV $ has a displacement generator ofa length at most $ \ell _T + \ell_V + 1 , $ which can be computed (over any field of constants) at thecost $ O ( \log n , n \log \log n )$.Furthermore, if $V$ is nonsingular, then the matrix $V^{-1} $ has a displacement rank  at most $ \ell_V + 2 $; this bound decreases to 2 if $V$ is a nonsingular Toeplitz matrix.}\vskip2mm{\bf Definition 2.3.} $\hbox{LIN} \cdot \hbox{SOLVE} (n) $ is theproblem of solving a nonsingular linear system $V \vec{x} = \vec{f}$ of$ n $ equations.  $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$ isa special case of $ \hbox{LIN} \cdot \hbox{SOLVE} (n) $, where thesystem has a Toeplitz-like coefficient matrix $V$ given with its displacement generator of a length $ \ell = O (1).$ The next definitions introduce an auxiliary problem related to$ \hbox{ LIN} \cdot \hbox{ SOLVE } \cdot \hbox{ T } (n).$ \vskip2mm{\bf Definition 2.4.}For a monic polynomial  $$C(x) = \sum_{i=0}^n c_{n-i}x^i = \prod_{k=1}^n (x-z_k ) \ , \ \ \ \ \c_0 = 1 \ , \eqno (2.1)$$ $s_j$ denotes the $j$-th power sum of its zeros, and $B(x)$denotes the reverse polynomial,$$s_j = \sum_{k=1}^n z_k^j \ , \ \ \ \ \ j=1,2,\ldots , \eqno (2.2)$$$$B(x) = \sum_{i=0}^n c_i x^i =x^n C(1/x) = \prod_{k=1}^n (1-xz_k ) \ .\eqno (2.3)$$ \vskip2mm{\bf Definition 2.5.}  a) $ \hbox{P} \cdot \hbox{SUMS} (n) $ denotes the problem of the computing of the power sums $ s_1 , \ldots , s_N $ of (2.2) for theinput set $ c_1 , \ldots , c_n $ of the coefficients of $B(x) $ and$C(x) $ of (2.1) and (2.3) and b) $ \hbox{I}\cdot \hbox{ P } \cdot \hbox{ SUMS} (n) $ denotes the problem of the inverse computation, of $ c_1 , \ldots ,c_n $ for the input $ s_1 , \ldots , s_N .$Here, $N=2^{\lceil \log m \rceil} <2m \leq 4n - 3,$ $m=n$ overthe fields of characteristic $0, m= n+\lfloor(n-1)/(p-1)\rfloor <2n$ over the fields of positive characteristic $p$. The solution of $\hbox{P} \cdot \hbox{SUMS} (n) $may rely onNewton's identities,$$s_k + \sum_{i=1}^{k-1} c_i s_{k-i} +kc_k = 0 \ , \ \ \ \k=1,\ldots ,n \ ,\eqno{(2.4)}$$$$s_k + \sum_{i=1}^n c_i s_{k-i} = 0 , \; \; k = n+1 , \ldots, N .\eqno{(2.5)}$$These identities form a triangular Toeplitz linear system of$ N $ equations in $ s_1 , s_2 , \ldots , s_N ; $ such a system canbe solved (over any field of constants) [BP94] at the cost bounded by$ O \left( ( \log n ) \log^{\ast} n , ( n \log \log n ) / \log ^{\ast} n\right) , $ since $ N < 4 n .$ In this paper we actually care only about $ \hbox{I} \cdot \hbox{P}\cdot \hbox{ SUMS} (n) $ because of the next result. \vskip2mm{\bf Fact 2.5} [P90/92]. $\hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T}(n) $ {\it for a linear system} $ V \vec{x} = \vec{f} $ {\it can bereduced (over any field of constants) to $\hbox{I} \cdot \hbox{P} \cdot\hbox{SUMS} (n) $ for $C(x) $ equal to the characteristic polynomial of$VW,$ for any fixed nonsingular Toeplitz-like matrix $W;$ the arithmetic costof this reduction is bounded by $O \left( ( \log n )^2 , ( n^2 \log \logn ) / \log n \right).$} {\bf Proof.} Fact 2.5 for $W = I$ is supported by the main algorithm of [P90/92].  The extension to any nonsingularToeplitz-like matrix $W$ follows from facts 2.3 and 2.4 since one mayfirst solve the linear system $(VW) \vec{y} = \vec{f} $ and then compute thevector $ \vec{x} =W \vec{y} , $ satisfying $V\vec{x} = \vec{f}$. \hfill $\Box$ The main algorithm of [P90/92] and its analysis in [P90/92] relying on fact 2.5 have constituted a major step towards the solution of the problem $\hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$. It remained to solve the problem $\hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n)$, and we will do this in sections 3--5. In these sections, we will rely on  Newton's identities of (2.4), (2.5) rewritten in a more compact form, as the polynomial identity   $$xB^{\prime} (x) =-B(x) \sum_{i=1}^{\infty} s_i x^i \ , \eqno (2.6)$$ which can be obtained from the next simple basic identity, $$xB^{\prime} (x)/B(x) =- \sum_{k=1}^nx z_k /(1-xz_k ) \ .$$\vskip2mm{\bf Definition 2.6.}  The $ n \times n $ matrix $ F = F ( \vec{c} ) =\left( \begin{array}{cccccc} & 0   & 1       &        & O    &   \\  &     & 0       & \ddots &      &   \\ &     &         & \ddots &      &   \\ & O   &       &        & 0    & 1 \\ &c_n  & c_{n-1} & \ldots & c_2  & c_1  \end{array} \right) $(having its first superdiagonal filled with ones, its last row with the vector $\vec{c} = (c_n , \ldots , c_1 ) , $ and all the other entries with zeros) is called the{\it Frobenius (companion) matrix} of the polynomial $ C(x) $ of(2.1). For a square matrix $T$, we write $C_T(x) = \det (xI-T)$ to denote the {\it characteristic polynomial} of $T$.   \vskip2mm{\bf Fact 2.6.} $C_F(x)= \det ( x I - F ) $,  {\it the characteristicpolynomial of } $ F = F ( \vec{c} ) $, has the coefficient vector $\vec{c}$ and, therefore, coincides with $C(x)$ of (2.1). \vskip2mm{\bf Fact 2.7.} {\it The matrix $ F$ has a displacement generator $(X,Y),\  \   X^T=\left( \begin{array}{ccccc}  1  & 0      & \ldots      & 0    &  0 \\   0  & 0      &  \ldots     & 0    & 1  %&c_n  & c_{n-1} & \ldots & c_2  & c_1  \end{array} \right)$, \\  $ Y^T=\left( \begin{array}{cccccc} 0    & 1       & 0      & \ldots      & 0   &  0 \\  c_n  & c_{n-1} & c_{n-2} &\ldots      & c_2  & c_1  \end{array} \right)$,of a length at most 2.} In section 5, we will use randomization based on thefollowing well-known fact, due to [DL78] (compare also [Z79], [S80], and [BP94], Lemma 1.5.1 on page 43). {\bf Fact 2.8.} {\it Let $P(x)=P(x_1 ,x_2 ,\ldots ,x_m )$ be a nonzero$m$-variate polynomial of a total degree $d$. Let $S$ be a finite subsetin the given field of constants or in its fixed extension field and let $X^* =\{ x_1^* ,\ldots ,x_m^* \}$be a point in $S^m$, where all $x_i$ are chosen at random in $S$,independently of each other, under theuniform probability distribution on $S$. Then${\rm Probability} \ \{P(X^* ) = 0 \} \leq d/|S|$,where $|S|$ denotes the cardinality of $S$.} {\bf Definition 2.7.} Hereafter, $S$ will denote a fixed (and sufficiently large) finite set, and all random values will be chosen from $S$ independently of each other and under the uniform probability distribution on $S$. Our randomized estimates for computational cost will not include the cost of generation of such random parameters. $|S|$ will denote the cardinality of $S$. The next definition is important for our algorithm of section 5 and is also used in section 12. {\bf Definition 2.8} [BGY80], [BP94]. For two nonnegative integers $k$ and $h$ and for a formal power series $Q(x)=\sum_{i=0}^{\infty} q_ix^i$ or, as an equivalent alternative, for a fixed $N \geq k+h$ and for a polynomial $Q_N(x)=\sum_{i=0}^N q_ix^i$, the $(k,h)$ entry of the {\it Pad\'{e} approximation table} (hereafter referred to as the {\it Pad\'{e} table}) is defined as a pair of polynomials $U(x)$ and $V(x) \neq 0$ of degrees at most $k$ and $h$, respectively, such that $U(x) -Q(x)V(x) = 0 \bmod x^{k+h+1}$.This entry {\it degenerates} if $U(x)$ and $V(x)$ have a nonconstant common divisor or if $\deg U(x)< k$ and $\deg V(x) <h$. {\bf Fact 2.9} [Gr72] (compare [BGY80], theorem 2). {\it Suppose that the polynomials $U(x)$ and $V(x)$ of definition 2.8 are relatively prime, that is, have no nonconstant common divisors. Then we have the 2 simultaneous bounds, $k > \deg U(x)$ and $h> \deg V(x)$, if and only if the associated $h \times h$ Toeplitz matrix $Q=(q_{i,j})$ is singular, where $q_{i,j}=q_{k-i+j}$, $i, j=0, 1, \ldots, h-1$, $q_\ell=0$ for $\ell<0$. Furthermore, $r={\rm rank}\  Q = h-\min ( k-\deg U(x), h-\deg V(x))$.} {\bf Corollary 2.1.} {\it Under the assumptions of fact 2.9, the $r\times r$ southeastern submatrix of the matrix $Q$ is nonsingular. }  {\bf Fact 2.10.} {\it The computation of the $(k, h)$ entry of the Pad\'{e} table for 2 given nonnegative integers $k$ and $h$ and for a formal power series $Q(x) = \sum_{i=0}^{\infty} q_ix^i$ can be reduced to computing the rank $r$ of an $h \times h$ Toeplitz matrix, solving a nonsingular Toeplitz linear system of $r$ equations, and multiplication of a $k \times k$ triangular Toeplitz matrix by a vector.  Furthermore, the rank computation can be omitted if it is known that the $(k,h)$ entry of the Pad\'{e} table does not degenerate.} Fact 2.10 is supported by Algorithm 2.5.1 from page 140 of [BP94], which relies on the proof of theorem 1 of [BGY80] and on our fact 2.9 and corollary 2.1.{\bf Definition 2.9.} The $(m+n) \times (m+n)$ matrix $R=[U,V]$, whose 2 blocks $U$ and $V$ are the $(m+n) \times m$ and $(m+n) \times n$ Toeplitz matrices, respectively, defined by their first rows $[u_n, 0, \ldots, 0]$ and $[v_m, 0, \ldots, 0]$, respectively, and by their first columns $[u_n, u_{n-1}, \ldots, u_0, 0, \ldots, 0]^T$ and $[v_m, v_{m-1}, \ldots, v_0, 0, \ldots, 0]^T$, respectively, is called the \\  {\it Sylvester or resultant matrix}, associated with the polynomials $u(x)=\sum_{i=0}^n u_i x^i$ and  $v(x)=\sum_{i=0}^m v_i x^i$. The $(m+n-2k) \times (m+n-2k)$ submatrix $R_k=[U_k,V_k]$ obtained by deleting the last $2k$ rows of $R$, $U$, and $V$ and the last $k$ columns of  each of $U$ and $V$ is called the $k$-th {\it subresultant matrix}, for $k=0,1, \ldots, \min \{ m, n \}$. (In particular, $R_0=R$.) $\det R$ is called the {\it resultant} of the 2 polynomials $u(x)$ and $v(x)$.By inspection, we immediately deduce the next result. {\bf Fact 2.11.} {\it The $k$-th subresultant matrix $R_k$ has a displacement generator of length $2$, given by the $(m+n-2k) \times 2$ matrices $X^{(k)}=\left( \begin{array}{lcr} u_n & \ldots & u_{2k-m+1} \\v_m-u_{n+m-k} & \ldots & v_{2k-n+1} -u_{k+1}\end{array} \right)^T$ and $Y^{(k)}=[\vec{e}^{(0)}, \vec{e}^{(m-k+1)}]$, where %$\vec{u}^h_g=[u_h, \ldots, u_g]^T$, %$\vec{v}^h_g=[v_h, \ldots, v_g]^T$, $u_k=v_k=0$ if $k<0$, $u_j=0$ if $j>n$, and the vector $\vec{e}^{(h)}$ has its $h$-th coordinate equal to 1 and other coordinates equal to 0.}{\bf Fact 2.12.} {\it Under definition 2.9, $u(x)$ and $v(x)$ have a nonconstant common divisor and/or $u_n=v_m=0$ if and only if $\det R = 0$.}{\bf Corollary 2.2.} {\it Under definitions 2.8 and 2.9, the $(k,h)$ entry of Pad\'{e} table degenerates if and only if the associated resultant vanishes, that is, $\det R = \det R(U, V)=0$, for $U(x)$ and $V(x)$ considered as polynomials of degrees $k$ and $h$, respectively.}  \noindent{\bf 3. From power sums to coefficients over the fields of characteristic $0$ or greater than $n$.}In this section we will solve the problem$ \hbox{I} \cdot \hbox{P} \cdot \hbox{ SUMS }(n) $ over the fields of constants that have characteristics $ p = 0 $ or $ p > n . $  In this case, the block substitution algorithm reduces the solution of the system (2.4) to multiplication of an $\lfloor n/2 \rfloor \times \lfloor n/2 \rfloor$Toeplitz matrix by a vector and to solving a pair of structured linear systems, similar to (2.4) but  of half-size [BP94]. By recursively applying blocksubstitution to the two latter systems and by using fact 2.3, we may compute $c_1 ,\ldots ,c_n$ at the overall cost $O(n\log n,(\log n)\log\log n)$. Actually, we will apply another algorithm, running faster (see theorem 3.1 at the end of this section). We have extracted the construction of this algorithm from the proof of lemma 13.1 on page 40 of [Sc82], where it was assumed that the input polynomial had all its zeros in the unit disc $\{ x: \  |x| < 1\}$ and where the presentation was directed towards estimating the bit-complexity of the solution. The algorithm is iterative and exploits the equation (2.6). The input is formed by the power sums of the zeros of $B(x)$, and the coefficients $c_0$ and $c_1$ of $B(x)$ are readily available. Each recursive step of the algorithm doubles the number of the available coefficients. Let us specify this algorithm. Write $$H_k (x) =B(x) \ {\rm mod} \ x^{k+1} = \sum_{i=0}^k c_i x^i \ ,\eqno (3.1)$$observe that $H_1 (x)=1-s_1 x$ is available, and recursively compute thecoefficients of $H_2 (x),H_4 (x),\ldots , H_{2^\nu} (x)$, $\nu = \lceil \log n\rceil$, by using the equations$$H_{2k} (x) =H_k (x) R_k (x) \ {\rm mod} \ x^{2k+1} \ , \eqno (3.2)$$where $R_k(x)$ are unknown auxiliary polynomials, and we easily deduce from (3.2) that $$R_k (x) = \sum_{j=0}^{2k} r_{k,j} x^j =1+x^{k+1} \sum_{i=0}^{k-1} r_{k,k+i+1} x^i\ ,\eqno (3.3)$$$k=1,2,4,\ldots ,2^{\nu-1}$. From (3.1)--(3.3), it is straightforward to deduce that$$R_k^{\prime} (x) = (R_k^{\prime} (x)/R_k (x)) \ {\rm mod}\ x^{2k+1} \ ,$$\begin{eqnarray*}H_{2k}^{\prime} (x)/H_{2k} (x) & = &(H_k^{\prime}(x) R_k (x) +H_k (x) R_k^{\prime} (x))/(H_k (x) R_k (x)) \\ & = & (H_k^{\prime} (x)/H_k (x)) +R_k^{\prime} (x) \ {\rm mod } \ x^{2k+1} \ .\end{eqnarray*} By combining the latter equation with (2.6), we obtain that$$-R_k^{\prime} (x)=H_k^{\prime} (x)/H_k (x) + \sum_{i=1}^{2k+1} s_i x^{i-1} \ {\rm mod} \ x^{2k+1} \ . \eqno (3.4)$$ Given $s_1 ,\ldots ,s_{2k+1} $ and $H_k (x)=B(x) \ {\rm mod} \ x^{k+1}$,we apply (3.4) in order to compute the coefficients of the polynomial$R_k^{\prime} (x) =\sum_{i=0}^{k-1}(k+i+1)r_{k,k+i+1}x^{k+i}$ [compare (3.3)].Then, we immediately recover the coefficients $r_{k,k+i+1}$ of$R_k (x)$, for $k+i+1\leq n$. (Under this bound, $k+i+1 \neq 0 \bmod p$for $p>n$,  so that $k+i+1$ has reciprocal in any field of characteristic0 or greater than $n$.) Finally, we compute $H_{2k} (x)$ by using (3.2). The entire computation amounts to computing polynomial sumsand products, as well as thereciprocals of polynomials with $x$-free terms equal to 1. More precisely, at the $k$-th recursive step, we perform $O(1)$ such operations modulo $x^{2k+1}$. We may implement them either  at the cost $O((\log k) \log^* k , (k \log \log k) / \log^* k)$ or, at our alternative choice, at the cost $O(\log k, ( \log \log k)^2 k)$ for $k=1, 2, 4, 8 \ldots, \lceil \log n \rceil -1$ (compare facts 2.1 and 2.2). By using Brent's principle, we slow down $ O(\log \log k)$ initial steps slightly, so as to be able to use fewer processors, and arrive at the following result. {\bf Theorem 3.1.} {\it The problem$\hbox{I} \cdot \hbox{P} \cdot \hbox{ SUMS } (n) $of definition 2.5 can besolved over any field {\bf F} of characteristic$p=0$ or $p>n$ at the cost bounded by$O((\log n)^2 ,(\log\log n)^2 n/\log n)$ or, alternatively, by$O((\log n)^2\log^* n,(n\log\log n)/((\log n)\log^* n))$. Moregenerally, the latter cost bound applies to the problem of theevaluation (over such a field {\bf F}) of the first $n$ coefficients ofany formal power series $B(x)$ satisfying the equations(2.6) modulo $x^{2n+1}$ and $B(0)=1$, where $s_1 ,\ldots ,s_{2n+1}$ are the input values. }%\vspace{.25in} \noindent{\bf 4. Extension to the computation over any field.} By extending the algorithm of the previous section to the computation over any field of constants that has a characteristic $p>0$, we may immediatelycompute at first the coefficients of $R_k^{\prime} (x)$ and then the coefficients $r_{k, k+i+1}$ of $R_k (x)$, except for ones where $p$ divides $k+i+1$, since $(r_{k,sp} x^{sp})^{\prime} =0$ mod $p$ for any integer $s$. Thus, for $0<p\leq n$, we need to extend the algorithm further in order to compute $B(x)$.Taking point of view of [KP92] and [KP94], we observe that the linear system of equations (2.4) has a matrix of rank $n - \lfloor n/p \rfloor$, and we need $\lfloor n/p \rfloor$ additional equations from (2.5) in order to obtain a nonsingular system. By associating the $k$-th equation of (2.4) with the power sum $s_k$, we may state the equivalent problem of finding a set $K$ of $n$ positive integers, $K=\{k_1, \ldots, k_n \}$, such that all the coefficients of the polynomial $B(x)$ can be rationally expressed through the set $s_{k_1}, \ldots, s_{k_n}$, called a {\it fundamental set} of the power sums. A. Sch\"{o}nhage in [S93] has shown (citing [V900], [K25], [K27], and [N27] as the preceding works) that the desired properties hold already for the set $K$ of the first $n$ positive integers not divisible by $p$.Then, based on this result, he has devised a parallel algorithm for the computation of the coefficients of the characteristic polynomial of a given general matrix over any field of constants. Of the 2 possible ways to devising such an algorithm, which both start with a given fundamental set of power sums of its zeros, he has chosen a way based (like in the preceding work [KP92]) on the reduction of the problem to solving a nonsingular linear system of $\lfloor n/p \rfloor$ equations. This way has lead him to an algorithm supporting estimates for the arithmetic parallel complexity of the solution that were strongly inferior to  ones based on the preceding algorithms of [Be84] and [Ch85]. In this section, we will recall his techniques, though  with more complete elaboration of the transition to the auxiliary power series $A(x)$; moreover, we will introduce randomization and will obtain improved estimates for parallel arithmetic complexity (see theorem 4.1, fact 4.1, and corollary 4.1). In the next section, we will demonstrate a  further crucial improvement, due to using an alternative approach, based on computing a Pad\'{e} approximation of $A(x)$. To make this  alternative approach work, we will relate it to the associated Frobenius (companion) matrix and to some other Toeplitz-like matrices and then will define a recursive process of computing Pad\'{e} approximation by exploiting the results of section 2 and some new techniques of regularization of Pad\'{e} approximation via randomization. This will lead us to parallel algorithms supporting the desired (improved) estimates for parallel arithmetic complexity. (Our presentation in these 2 sections have some similarity with the presentation in the paper [P93a], partially reproduced on pp. 373--377 of [BP94], where similar techniques have been studied   for computations with general matrices.)We first rewrite $B(x)$ as follows:$$B(x) = \sum_{k=0}^{p-1} B_k (x^p )x^k \ , \eqno (4.1)$$where$$B_k (y) = \sum_{i=0}^{d_k} c_{k+ip} y^i \ , \ \ \ \ \k=0,1, \ldots ,p-1 \ , \eqno (4.2)$$$$d_0 =d=\lfloor n/p\rfloor \ , \ \ \ \d_k \geq d_{k+1} \geq d-1 \ ,  \ \ \ \k=0,1, \ldots ,p-2 \ . \eqno (4.3)$$Then we observe that $B_0 (0) =B(0)=1$ and define $p-1$ rationalfunctions,$$A_k (y) =B_k (y)/B_0 (y) \ , \ \ \ \ \ k=1,\ldots ,p-1 \ , \eqno (4.4)$$as well as the next rational function, which we also represent as aformal power series:$$A(x)=B(x)/B_0 (x^p )=\sum_{i=0}^{\infty} a_i x^i \ . \eqno (4.5)$$Combining (4.1), (4.2), (4.4), and (4.5), we obtain that$$A(x) =1+ \sum_{k=1}^{p-1} x^k A_k (x^p ) \ , \eqno (4.6)$$which implies that$$a_0 = 1 \ , \ \ \ \ \ a_{sp} = 0 \ , \ \ \ \ \ s=1,2,\ldots , \eqno (4.7)$$$$a_k =c_k \ , \ \ \ \ \ {\rm for} \ k=1,2,\ldots ,p-1 \ . \eqno (4.8)$$(4.5) and (4.6) also imply representation of $A_k(y)$ as formal power series for all $k$. Now recall (4.2) for $k=0$ and observe that $B_0^{\prime} (x^p )=0$mod $p$. Therefore, (4.5) implies that $A^{\prime} (x) =B^{\prime} (x)/B_0 (x^p )$, $x  A^{\prime} (x)/A(x)=xB^{\prime} (x)/B(x)$. Combine the latter equation with (2.6) and obtain that$$xA^{\prime} (x)=-A(x) \sum_{k=1}^{\infty} s_k x^k \ . \eqno (4.9)$$ We will now prove the following result. {\bf Theorem 4.1.} {\it Given the power sums $s_k$ of the zeros of apolynomial $C(x)$, $k=1, 2, \ldots, 2n$ [see (2.1), (2.2)], one may compute (over any field ofconstants of a positive characteristic $p$) the coefficients $a_1 ,a_2 ,\ldots ,a_m, \ m=n+\lfloor (n-1)/(p-1)\rfloor < 2n$,of the formal power series $A(x)$ satisfying (4.5), (4.7) [for$B(x)$ of (2.3) and $B_0 (x^p)$ of (4.1)] at the parallel arithmetic cost bounded by $O((\log n)^2, (n(\log\log n)^2 )/\log n)$ or, alternatively,$O((\log n)^2\log^* n,(n\log\log n)/((\log^* n)\log n))$, which are thesame bounds as in theorem~3.1.} {\bf Proof.} The computation starts with evaluating $a_i =c_i$ for$i<p$ [see (4.8)], by means of the algorithm of the previous section.Then we set $a_p =0$, $k=p$, and again apply the algorithm of the previoussection, this time, however, replacing $c_i$ by $a_i$ in (3.1), $B(x)$ by$A(x)$ and $B^{\prime} (x)$ by $A^{\prime} (x)$ throughout, so that, inparticular, (2.6) is replaced by (4.9), and $H_k(x)$ is set to equal $\displaystyle{\sum_{i=0}^k}a_ix^i$. The resulting algorithm works asbefore, except that $r_{p,2p}$ cannot be recovered from $R_p^{\prime} (x)$but is now computed from the equations$$-\sum_{i=0}^{p-1} a_i r_{p,2p-i} = a_{2p} \ ,  \eqno (4.11)$$implied by (3.2) and (3.3) for $k=p$. Substituting (4.7), we deduce that$$r_{p,2p} =-\sum_{i=1}^{p-1} a_i r_{p,2p-i} \ .$$Then we immediately compute $r_{p, 2p}$. The same process is recursivelyrepeated for $k=2p,4p,8p,\ldots ,$ until all the $a_i$, $i=1,2,\ldots ,m$,have been computed. Each recursive step begins with the computation of thecoefficients $r_{k,sp}$ of $R_k (x)$ for $k=qp$, $s=q+1,\ldots ,2q$. Thiscomputation relies on the equations$$-\sum_{i=0}^{sp-1} a_i r_{qp,sp-i} =a_{sp} \ , \ \ \ \ \s=q+1,\ldots ,2q \ ,$$which extend (4.11) from the case $s=2,\ q=1$. Substitute (4.7)and obtain the equations$$r_{qp,sp} =- \sum_{i=1}^{sp-1} a_i r_{qp,sp-i} \ , \ \ \ \ \s=q+1, \ldots ,2q \ .$$Now replace $r_{qp, sp-i}$ by $r_{qp, sp-i}^*$ on the right sides, where$r_{qp, sp-i}^* = 0$ if $p$ divides $i$, $r_{qp, sp-i}^* =r_{qp, sp-i}$ otherwise.Since $a_i=0$ if $p$ divides $i$ and if $i> 0$, the right sides do not change,and therefore,$$r_{qp, sp} = - \sum_{i=1}^{sp-1} a_i r_{qp, sp-i}^* \ , \ \ \ \ \s=q+1, \ldots ,2q \ . \eqno (4.12)$$Here, $a_i$ and $r_{qp, sp-i}^*$ are known for all $i$, and we obtain all the$r_{qp, sp}$ from (4.12) as the coefficients of a polynomial product, at the cost$O(\log q,q\log\log q)$, (compare fact 2.1). Having computedthe $r_{qp, sp}$ for all $s$, we then recall (3.2) and obtain thecoefficients $a_{qp+1} ,\ldots ,a_{2qp}$ of $H_{2qp} (x)$, at the cost withinthe same bound, due to fact 2.1. This completes the recursive step ofcomputing the coefficients $a_{k+1},a_{k+2}, \ldots ,a_{2k}$ for $k=qp$. Thecomputational cost of this step is bounded according to theestimates of section 3 for the cost of the transitionfrom $H_k(x)$ to $H_{2k}(x)$ for $k=qp$. Therefore, in$\lceil \log n\rceil +2$ recursive steps, we compute $a_1,\ldots, a_{2n}$at the cost bounded as we claimed in theorem 4.1. \hfill $\Box$ Having computed the first $m$ coefficients $a_1 ,\ldots, a_m$, for$m=n+ \lfloor (n-1)/(p-1) \rfloor < 2n$, we then obtain the coefficientsof $B(x)$ as follows: \vskip2mm{\bf Algorithm 4.1.} {\bf Input:} $a_1, \ldots, a_m$. {\bf Output:} $c_0, c_1, \ldots, c_{n-1}$.{\bf Computations:}\begin{enumerate}\item Compute the coefficients $c_p ,c_{2p} ,\ldots ,c_{dp}$ of$B_0 (y)$ [see (4.1)-(4.3)];\item Compute the remaining coefficients of $B(x)$.\end{enumerate} Due to (4.4)and (4.6), stage 2 immediately reduces to $p-1$ multiplications of pairs ofpolynomials modulo $x^{d+1}$, $d=\lfloor n/p \rfloor$,[see (4.1)-(4.3)], andmay, therefore, be performed at the overall cost $O(\log (n/p),n \log\log(n/p))$, as this follows from fact 2.1. Stage 1 can be reduced to solving a linear system of $d= \lfloor n/p \rfloor $ equations in $c_p, c_{2p}, \ldots, c_{dp}$. The $d \times d$ matrix $Q$  of this system is filled with the coefficients $a_i$ of $A(x)$ for $i=0, 1, \ldots, m$, where $m = m_{n,p}=n+\lfloor (n-1)/(p-1) \rfloor < 2n $, as in theorem 4.1. {\bf Fact 4.1.} {\it $\det Q \neq 0$ with a probability at least  $1 - (1+\lfloor m/p \rfloor )m /|S|$, $m<2n$, provided that  the coefficients of $B(x)$ are chosen from a finite set $S$ (cf. definition 2.7), independently of each other and under the uniform probability distribution on $S$.}  {\bf Proof.} In [Sc93], A. Sch\"{o}nhage has proved nonsingularity of $Q$ in the case of generic polynomial $B(x)$. [He gave 2 proofs. In one of them, he showed that the diagonal entries of $Q$, considered as indeterminates, made contribution to $\det Q$ that was unmatched by all other terms. In another proof, he specified that $Q$ turned into the $m \times m$ identity matrix (having its determinant equal to 1), under the choice of $C(x) = 1+\sum_{i=1}^{p-1} x^{i+r(i)p}$, for some integers $r(i)$ satisfying $ -1 \leq \lfloor (n-1)/(p-1) \rfloor - \lfloor n/p \rfloor-r(i) \leq 1$, $i=1, \ldots, p-1$.] To complete the proof of fact 4.1, it remains to show that $\det Q$ is a polynomial of a degree at most $(1+ \lfloor m/p \rfloor)m $ in the coefficients $c_1, c_2, \ldots, c_n$ of $C(x)$ and $B(x)$ (and then, fact 4.1 will follow from fact 2.8). Let us do this. Recall (4.5) and obtain that $A(x) \bmod x^{m+1} = B(x)\sum_{i=0}^{\infty} (1-B_0(x^p))^i \bmod x^{m+1}$. Since $x^p$ divides $1-B_0(x^p)$, we have $(1-B_0(x^p))^i \bmod x^{m+1} =0$ for $i>m/p$. Therefore, the coefficients $a_1, a_2, \ldots, a_m$ of $A(x) \bmod x^{m+1}$ are polynomials in $c_1, c_2, \ldots, c_n$ of degrees at most $1+\lfloor m/p \rfloor$. Since $Q$ is an $m \times m$ matrix, $\det Q$ is a polynomial in $c_1, c_2, \ldots, c_n$ of a degree at most $(1+\lfloor m/p \rfloor)m$. \hfill $\Box$ Let us now assume that the coefficients $c_1, c_2, \ldots, c_n$ of $C(x)$ and $B(x)$ have been randomly chosen (from a fixed finite set $S$ of definition 2.7) but remain unknown, whereas the power sums $s_1, s_2, \ldots, s_m$ are known.  Such an assumption can be implicitly motivated by fact 4.1, but our actual motivation is more direct. It  relies on the reduction of our main problem $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$ for an arbitrary nonsingular Toeplitz-like input matrix $V$ to the same problem but for a random Toeplitz-like input matrix $VW$ and, consequently, to the associated problem $\hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n)$, corresponding to a random characteristic polynomial $C(x)$. We will exploit such a reduction in the next section.  By summarizing our analysis of stages a) and b) of algorithm 4.1, we arrive at the following result. {\bf Corollary 4.1.} {\it For a pair of natural $n$ and $p$, for$m=n+\lfloor (n-1)/(p-1)\rfloor <2n$, and for a given set of the values $s_1, s_2, \ldots, s_m$ of the power sums of the zeros of a random polynomial $C(x)$ of (3.1),with its unknown coefficients (chosen from a fixed set $S$ of a cardinality $|S|$, independently of each other and under the uniform probability distribution on $S$),the problem $ \hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n) $of definition 2.5 can be solved (over any field ofconstants of a positive characteristic $p$),with a probability at least $1 - (1+\lfloor m/p \rfloor) m/|S|$,at the asymptotic cost boundedby the cost bound of theorem 4.1 augmented by the complexity of solving the problem $\hbox{LIN} \cdot \hbox{SOLVE}(\lfloor n/p \rfloor)$.}\vspace{.25in} \noindent{\bf 5. Improved algorithms and complexity bounds.} In this section, some new techniques will enable us to  extend corollary 4.1 by devisinga randomized algorithm for theproblem $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n) $ over anyfield of constants of a characteristic $p>0$.Moreover, we will simplify the computation of thecoefficients of $B_0(y)$ and will arrive at the improved complexityestimates. Due to (4.4), the coefficients $c_p ,c_{2p} ,\ldots ,c_{dp}$ of $B_0 (y)$can be obtained from the $(d_1 ,d)$ entry of thePad\'{e} approximation table for the power series $A_1 (y)$ (compare definition 2.8),where $A_1(y)=A_1 (x^p)$ is defined by $A(x)$, due to (4.5) and (4.6), and $d_1 \leq d=\lfloor n/p \rfloor$, due to (4.3). Unless degeneration occurs,that is, unless $c_n = c_{n-1} = 0$ and/or $ B_1 (y) $ and$ B_0 (y) $ have a nonconstant common divisor, such a$(d_1 ,d)$ entry of the table is filled with the pair $(B_1 (y) ,B_0 (y))$ and can be computedby means of the known reduction of the problemto solving a nonsingular Toeplitz linear system of $d$ equations(see fact 2.10). We will first apply randomization to make degeneration unlikely.  We will call this technique {\it regularization of Pad\'{e} approximation via randomization}. \vskip2mm{\bf Fact 5.1.}  Let $V $ {\it denote a nonsingular $ n \times n $Toeplitz or Toeplitz-like matrix given with its displacement generator of length $ \ell$.Let $X$ and $Y$ be a pair of $ n \times ( \ell + 5 ) $ matriceswith random entries (chosen from a fixed finite set $S$of cardinality $ | S | , $ independently of each other and under the uniform probability distribution on $S$).Let $W$ be the Toeplitz-like matrix defined by its displacement generator$X,Y.$  Let $C (x) $ denote the characteristic polynomial of the matrix$ F = VW , $ and let $ B(x), B_0 (y)$ and $ B_1 (y) $ be the associatedpolynomials of (2.3) and (4.2) for an integer $ p> 1 .$  Then, \\a) the probability of degeneration of the pair $ B_1 (y) ,B_0 (y) $ is at most $ 4 d^2n / | S| , \;d =\lfloor n/p \rfloor$, and \\b) the probability that the matrix $W$  is singular is at most $2n / |S|$.} \vskip2mm{\bf Proof.}  Due to fact 2.9, the degeneration implies singularity of the $d \times d$ Toeplitz matrix $A^{(1)}=(a_{i,j}^{(1)})$, $ a_{i,j}^{(1)}=a_{d_1-i+j}^{(1)}$; $i,j=0, 1, \ldots, d-1$, $a_h^{(1)}=0$ for $h<0$, associated with the formal power series $A_1(y)= B_1 (y)/B_0 (y) = \sum_{i=0}^{\infty} a_i^{(1)} y^i$, where $a_i^{(1)}=a_{ip+1}$, $i=0, 1, \ldots$ [cf. (4.5), (4.6)]. The determinant of this matrix is a polynomial  of a degree at most $ d $ in the coefficients of $A_1 (y) $ and, therefore, is a polynomial of a degree at most $(d+d_1)d \leq  2d^2$ in the coefficients of $C(x).$  These coefficients are inturn polynomials of degrees at most $n$ in the entries of $VW$ and, therefore, also in the entries of $W.$  The latter entries are bilinear in the randomentries of $X$ and $Y.$  Therefore, the degeneration of the pair$B_1 (y) , B_0 (y) $ implies the vanishing of some polynomial of a degree at most $4d^2n $ in the entries of $X$ and $Y.$  Due to fact 2.8, the part a) of fact 5.1 will follow as soon as we show that the associated pair $B_1 (y) ,B_0 (y) $ does not degenerate for some special choice of $X$ and $Y$. Now, choose any nondegenerating pair of polynomials of degrees $d_1 $ and $ d,$ denote them $ B_1 (y) $ and $B_0 (y),$ extend them to some pair of polynomials$B(x)$ and $C(x)$ satisfying (2.3), (4.1)-(4.3), and define the associated Frobenius (companion) matrix $ F = F (\vec{c}).$  Finally, set $W=V^{-1} F$, so that $VW = F (\vec{c}),$ and the associated pair$ B_1 (y),B_0(y)$ does not degenerate.  Due to facts 2.4 and 2.7, the matrix $W$ has a displacement generator of a length at most $ \ell + 5 , $ which we will choose as the desired pair $ X , Y .$ This proves part a) of fact 5.1. Part b) follows from fact 2.8 since the entries of $W$ are quadratic polynomials in the entries of $X$ and $Y$ (see definition 2.2) and, therefore, $\det W$ is a polynomial of a degree at most $2n$ in these entries.\hfill $\Box$ Now, we are ready to define the following algorithm for $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$ over any field of constants of any characteristic $p.$ The algorithm will rely on a {\it bipartite recursive process of computing Pad\'{e} approximations and solving Toeplitz linear systems}, with recursive decrease of the input size by a factor at least $p$ in each step and with recursive regularization of  Pad\'{e} approximations via randomization. \vskip2mm{\bf Algorithm 5.1.}{\bf Input:} a field ${\bf F}$ of a characteristic $p >0$, a displacement generator of a length $\ell =O(1)$ for an $n \times n$ Toeplitz-like matrix $V$, an $n$-dimensional vector $\vec{f}$, and a set $S$ of a cardinality $|S|$ consisting of elements of the given field ${\bf F}$ or of its algebraic extension ${\bf E=E(F)}$.  {\bf Output:} either FAILURE, with a probability at most $2(2n+p)/((p-1)|S|)$, or a vector $\vec{x} = V^{-1} \vec{f}$, with a probability at least $1 - 2(2n+p)/((p-1)|S|)$.  {\bf Computations.}\begin{enumerate} \item{} %no 1Compute 2 random matrices $ X$ and $Y$ (with entries from $S$) satisfying the assumptions of fact 5.1. \item{} %no 2Compute a displacement generator of a length at most $ 2 \ell + 6 $ for the  matrix $VW$, where the matrix $W$ is defined by its displacement generator $ X , Y$ (see fact 2.4).\item{} %no 3Reduce the solution of the linear system $VW \vec{y} = \vec{f} $ to theproblem $ \hbox{I} \cdot \hbox{P} \cdot \hbox{SUMS} (n) $ for$ C (x) = \det (x I - VW ) , $ (by applying fact 2.5 and the main algorithm of [P90/92]).\item{} %no 4 If $ p = 0 $ or $ p > n , $ solve the latter problem $ \hbox{I} \cdot\hbox{P} \cdot \hbox{SUMS} (n) $ by means of the algorithm of section 3(see theorem 3.1) and output the solution to $ \hbox{LIN} \cdot\hbox{SOLVE} \cdot \hbox{T} (n).$  Otherwise, extend this algorithm to thecomputation of the coefficients $ a_1 , \ldots , a_m $ of (4.5), (seetheorem 4.1), and then go to stage 5.\item{} %no 5Apply algorithm 4.1 but at its stage 1 compute the coefficients of the polynomial $ B_0 (y) $ from the$ (d_1 , d ) $ entry of the Pad\'e approximation table for the formal power series $A_1 (y) $ defined by $A(x) $ [see (4.4), (4.6), and definition 2.8].  The entry $ (d_1 , d) $ of the Pad\'e table is computed (based on fact 2.10 and on the observation that the pair $B_1(x)$, $B_0(x)$ does not degenerate) via solving a nonsingular Toeplitz linear system of at most $ d = \lfloor n/p \rfloor $ equations.For this problem, $ \hbox{ LIN } \cdot \hbox{SOLVE} \cdot \hbox{T} (d)$,  algorithm 5.1 is recursively invoked until (in at most $g=1 + \lfloor (\log n)/ \log p \rfloor$ its recursive applications) the original problem is reduced to solvinga nonsingular Toeplitz linear system of at most $ p-1 $ equations;then a singleapplication of stage 4 of this algorithm defines a vector $\vec{x}$, which is a candidate for being the solution vector $V^{-1}\vec{f}$ to the original problem. (Bipartite transition between the Pad\'{e} entry computation and $ \hbox{T} \cdot \hbox{LIN} \cdot \hbox{SOLVE} (n) $ at every recursive step motivates our use of the name of bipartite recursion.)\item %no 6Compute the vector $V\vec{x}$ and compare it with $\vec{f}$. If $V\vec{x}= \vec{f}$, output $\vec{x}$. Otherwise output FAILURE.\end{enumerate} Let us apply and extend the preceding analysis in order to show {\it correctness } of the algorithm. Clearly, the algorithm is of Las Vegas type: it may output FAILURE but otherwise produces correct answer. Failure may only occur in the cases where one or more of the auxiliary matrices $W$  are singular and/or one or more of the auxiliary Pad\'{e} approximation problems degenerate, due to an unfortunate choice of random parameters. Let us estimate the failure probability and the number of random parameters involved.  {\bf Proposition 5.1.} {\it Algorithm 5.1 involves less than $(2\ell + 10 + 14/(p-1))n \leq (2 \ell +24)n$ random parameters.  }{\bf Proof.} Recall that at most $g=1 + \lfloor (\log n)/\log p \rfloor$ random Toeplitz-like matrices are involved in computations by algorithm 5.1; furthermore, due to fact 5.1, their displacement generators have lengths at most $\ell +5$ at the first recursive step and at most 7 at all the subsequent steps (since at the latter steps, we deal with Toeplitz linear systems arising from Pad\'{e} approximation problems and since for such systems we may replace $\ell$ by 2 in the estimate of fact 5.1). Therefore, the displacement generator is defined by $2(\ell +5)n$ parameters at the 0-th recursive step and by $14n_k \leq 14n/p^k$ parameters at the $k$-th recursive step, where $p \geq 2$, $n_k = \lfloor n_{k-1}/p\rfloor$, $n_0=n\ , \ \  k=1,2, \ldots$. This implies that a total of $2(\ell +5)n + 14\sum_{k=1}^{g-1} n_k \leq 2(\ell +5)n + 14n\sum_{k=1}^{g-1} p^{-k} < (2 \ell +10 +14/(p-1))n \leq (2\ell +24)n$ parameters are involved in algorithm 5.1. \hfill $\Box$ {\bf Proposition 5.2.} {\it Algorithm 5.1 outputs either FAILURE, with a probability at most $ \frac{2np}{|S|}(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)})$, or the correct solution.} {\bf Proof.} We recall that algorithm 5.1 fails only if, for some $k$, at the $k$-th recursive step of its computation, the auxiliary matrix $W$ is singular and/or the associated Pad\'{e} approximation problem degenerates. Due to fact 5.1, such 2 events may occur with probabilities bounded by $2n_k /|S|$ and $4n_k^3/(p^2|S|)$, respectively. Therefore, while we perform algorithm 5.1, we will encounter no singularity with a probability at least $$ \prod_{k=0}^{g-1} (1-\frac{2n_k}{|S|}) \geq 1 - \frac{2}{|S|}\sum_{k=0}^{g-1} n_k \geq 1- \frac{2n}{|S|}\sum_{k=0}^{g-1} p^{-k} > 1- \frac{2np}{(p-1)|S|} $$ and will encounter no degeneracy with a probability at least  $$ \prod_{k=0}^{g-1} (1-\frac{4n_k^3}{p^2|S|}) \geq 1 - \frac{4}{p^2|S|}\sum_{k=0}^{g-1} n^3_k \geq 1- \frac{4n^3}{p^2|S|}\sum_{k=0}^{g-1} p^{-3k} > 1- \frac{4n^3p}{(p^3-1)|S|} \ . $$ Therefore, algorithm 5.1 outputs the correct solution with a probability at least  $$(1- \frac{2np}{(p-1)|S|})(1- \frac{4n^3p}{(p^3-1)|S|}) > 1- \frac{2np}{(p-1)|S|} -\frac{4n^3p}{(p^3-1)|S|} = 1-\frac{2np}{|S|}(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)})\ .$$ It follows that the failure probability is at most $\frac{2np}{|S|}(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)})$. \hfill $\Box$ This completes the {\it correctness} proof for algorithm 5.1, and proposition 5.1 bounds the number of random parameters involved in this algorithm. By using our analysis in this and previous sections (in particular, cf. fact 2.10), we easily estimate the {\it computational cost} of performing algorithm 5.1, in the case where $|S|\leq |{\bf F}|$ [see (5.3) below]. We also need to consider the opposite case, where $|S|> |{\bf F}|$. Indeed, due to proposition 5.2, it suffices to choose $$ |S|= 2^{\beta} \lambda, \ \  \lambda=\lambda(n, p)= 2np(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)}) \ , \eqno (5.1) $$ in order to bound the probability of failure of algorithm 5.1 by $2^{-\beta}$, for any fixed $\beta$. The value $|S|$ of (5.1), however, may exceed  the cardinality $|{\bf F}|$ of the field ${\bf F}$ of constants (where $|{\bf F}| \geq p$ but can indeedbe as small as $p$). Then, we will choose the set $S$ in an algebraic extension field ${\bf E}$ of polynomials over ${\bf F}$ reduced modulo a fixed (irreducible) polynomial of degree $$ D=\lceil (\log |S|)/\log |{\bf F}| \rceil \eqno (5.2) $$(compare our remark 5.1). Such an extension field ${\bf E}$ has a cardinality $ |{\bf E}| \geq |S|$. We note that for $\beta= O(\log n)$, (5.1) implies that $\log |S|=O(\log n)$, $$D=O((\log n)/\log |{\bf F}|) \ . %\eqno (5.2a)$$Any $(k-1)$-variate polynomial over ${\bf E}$ can be represented as a $k$-variate polynomial over ${\bf F}$ whose degree in the new variable (defining the extension field ${\bf E}$ over ${\bf F}$) is at most $D-1$. On the other hand, inspection of algorithm 5.1 and of the main algorithm of [P90/92] (supporting stage 3 of algorithm 5.1) shows that all their operations are ultimately reduced to computing polynomial sums, products, and reciprocals modulo some powers of $x$. We combine the latter observation with (5.1), (5.2), remark 2.1, propositions 5.1, 5.2, and various results involved in the description of algorithm 5.1, and arrive at the following theorem. {\bf Theorem 5.1.} {\it Let ${\bf F}$ be a field of a characteristic $p$, let $\beta$ be a positive scalar, $\lambda = 2np(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)}) \leq 4n (1 + \frac{2n^2}{7})$ for $p\geq 2$ [cf. (5.1)], $N=2^{\beta}\lambda$, $D=\lceil (\log N)/ \log |{\bf F}| \rceil$, and $S$ be a fixed finite set of cardinality $|S|=N$ and with elements either from the field ${\bf F}$ if ${\bf F}$ has a cardinality $|{\bf F}|$ at least $N$ or from a field ${\bf E}$ being an algebraic extension of ${\bf F}$ of degree $D$. Then, algorithm 5.1 solves the problem  $ \hbox{LIN} \cdot\hbox{SOLVE} \cdot \hbox{T} (n)$ at a deterministic cost bounded by $C_n =O((\log n)^2, (n^2 \log \log n)/\log n)$, if $p=0$ or $p>n$ [compare (1.2)]. If $0< p \leq n$, algorithm 5.1 solves the same computational problem at a randomized cost bounded by }$$C_{n,p} = O((\log n)^2 \gamma (n,p),(n^2 \log\log n)/(\gamma (n,p)\log n ))   \ \   for \  \  |{\bf F}| \geq N  \eqno (5.3)$$  {\it and by} $$C_{n,p, {\bf F}} = O((\log n)^2 \gamma (n,p),(n^2 \log\log n)(\log N)/(\gamma (n,p)(\log n)\log |{\bf F}| ))  \ \  for \  \  |{\bf F}| < N  \eqno (5.4)$${\it [which turns into} $$C_{n,p, {\bf F}} = O((\log n)^2 \gamma (n,p),(n^2 \log\log n)/(\gamma (n,p)\log |{\bf F}| ))  \  \   if  \  \beta=O(\log n) ]  \ , $${\it where }$$ \gamma (n,p) = \lceil (\log n)/\log p \rceil \ . \eqno (5.5)$$ {\it  %[compare (1.3)]. In the cases where $0<p\leq n$, the computation involves at most $(2 \ell +10 +\frac{14}{p-1})n \leq (2 \ell +24)n$ random parameters (for $p \geq 2$), chosen from the set $S$, independently of each other and under the uniform probability distribution on $S$, where $\ell$ denotes the length of the displacement generator  that defines the Toeplitz-like input matrix. Furthermore, in the case where $0< p \leq n$, the output is either FAILURE with a probability at most $\frac{2np}{N}(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)})\leq 2^{-\beta}$ or a correct solution, whose correctness test has a cost $O(\log n, n \log \log n)$. The randomized cost estimates include neither the cost of computing the polynomial $E(x)$ that defines the algebraic extension field ${\bf E}$ (see remark 5.1 below) nor the cost of generating random parameters.}  {\bf Remark 5.1.} In theorem 5.1 and similarly in the theorems of the next sections, the degree $D$ of the auxiliary algebraic extension field is defined by (5.2) and is $O(\beta + \log n)$. If $\beta = O(\log n)$ and  ${\bf F} ={\bf F}_p$ is a finite field with $p$ elements, then the algorithm of [Sh94] computes a basic irreducible polynomial defining an algebraic extension field of degree $D$ over ${\bf F}$ at the cost $O((\log n)\log(np), 1)$, which for us turns into $O((\log n)^2,1)$ since $p\leq n$ in applications in this paper. {\bf Remark 5.2.} In algorithm 5.1, we may modify stages 2 and 3, so as to  choose random matrices $X$ and $Y^T$ of size $n \times (\ell +2)$ and a random $n \times n$ Frobenius matrix $F$ of definition 2.6, at stage 2, and to replace $VW$ by $VWF$, at stage 3. Then, only $(2\ell +19)n$ random parameters are needed in the algorithm, but our upper bounds on the probability of encountering singularity and/or degeneracy increase a little, so as to reach $\frac{2np+1}{(p-1)|S|}$ and $\frac{6n^3p}{(p^3-1)|S|}$, respectively. Similar modifications of randomization are possible in the  algorithms of the next sections. {\bf Remark 5.3.} Algorithm 5.1 for $W=I$ can serve as a heuristic algorithm for $\hbox{I}\cdot\hbox{P}\cdot\hbox{SUMS}(n)$, which may only fail in the case of degeneracy of a Pad\'{e} approximation  computed at stage 5 (also compare remark 7.3 in section 7).{\bf Remark 5.4.} The precise estimates for the number of random parameters and for the failure probability in theorem 5.1 and in the theorems of the next sections are quite complicated. The reader, however, may greatly simplify these estimates by writing them in a less specific form as follows: choosing $n^{O(1)}$ random parameters from a set $S$ of a cardinality $|S|= n^{O(1)}$ suffices in order to ensure simultaneously the potential work bound $O((n^2 \log  \log n)(\log n)^2/\log \min \{ |{\bf F}|, n \})$,   the time bound $O((\log n)^g (\gamma(n,p)^h)$, and the upper bound $2^{-\beta}$ for any positive $\beta = O(\log n)$ on the failure probability, where $|{\bf F}|$ is the cardinality of the given field of constants, $2\leq g \leq 3$, $0\leq h \leq 2$, and $\gamma(n,p)$ is defined by (5.5). The estimates in this form are easier to comprehend, in many cases they are sufficient, and their derivation is less tedious. \vskip2mm\noindent{\bf 6. Extension to computation of the inverses and the determinants of Toeplitz-like matrices.}Let an $n \times n$ Toeplitz-like matrix $V$ be given with its short displacement generator of a length $\ell = O(1)$.  Then the algorithms of [P90/92] enable us to compute (over a field of constants of a characteristic $p$) the characteristic polynomial $C_V(x)$, the determinant $\det V$, and if $\det V \neq 0$ [that is, if $V$ is nonsingular], then also $V^{-1}$, at the overall deterministic  cost $C_n$ of (1.2), provided that $p=0$ or $p>n$. Now, let $0< p \leq n$.  Further inspection of algorithm 5.1 immediately enables us to extend the estimates of theorem 5.1 to the problem of the computation (over any field ${\bf F}$) of the coefficients of the characteristic polynomials $C_{VW}(x)$ and  $C_W(x)$ of $n \times n$ matrices $VW$ and $V$, where $V$ and $W$ are $n \times n$ Toeplitz-like matrices, $V$ is nonsingular, a displacement generators of a length $\ell$ for $V$ is given as a part of the input, and a displacement generator of length $\ell +5$ is randomy chosen for $W$. The $x$-free terms of $C_{VW}(x)$ and $C_W(x)$ are equal to $(-1)^n\det (VW) = (-1)^n (\det V) \det W$ and $(-1)^n \det W$, respectively. With a probability at least $1-n/|S|$, we have $\det W \neq 0$ (due to fact 2.8), and then we obtain $$\det V = C_{VW}(0)/C_W(0) \ .   $$ The computational cost bounds of theorem 5.1 apply here too (though with the double number of random parameters and the double failure probability since we apply theorem 5.1 twice, to $VW$ and $V$). Furthermore, given a nonsingular $n \times n$ Toeplitz-like matrix $T$ and its characteristic polynomial $C_T(x)$, we may apply an algorithm from [P90/92] (also available in [BP94]) in order to compute  the matrix $T^{-1}$ at the cost bounded by $C_n$ of (1.2). By applying this result to $T=VW$, we obtain $(VW)^{-1}$ and $V^{-1}=W(VW)^{-1}$ at the overall randomized cost bounded according to theorem 5.1. The presented randomization is of Monte Carlo type. To yield Las Vegas randomization for $V^{-1}$, we multiply, at the cost $O(\log n, n^2\log \log n)$, the computed matrix $V^{-1}$ by the input matrix $V$. If the product gives us the identity matrix, then the matrix $V^{-1}$ has been computed correctly, and we output it. Otherwise, we output FAILURE.  Clearly, the algorithm always outputs FAILURE if $V$ is singular, and our randomized computation of $V^{-1}$ is of  Las Vegas type, indeed.  On the other hand, we recall that randomization was only needed in the transition from the power sums $s_i$ to a characteristic polynomial. Therefore, we may alternatively verify correctness of the output matrix $V^{-1}$ (as well as of the output value $\det V$) by testing correctness of such transitions to   the characteristic polynomials $C_{VW}(x)$ and $C_W(x)$. The latter tests are reduced to solving triangular Toeplitz linear systems of Newton's identities in the power sums $s_i$ [cf. (2.4), (2.5)], which can be done at the cost$O((\log n) \log^*n, (n\log \log n)/\log^* n)$ [BP93]. Furthermore, the same correctness tests can be applied to certify nonsingularity of $V$. [$V$ is certified to be nonsingular unless FAILURE is output. In the latter case, we may conclude (with a low error probability) that $V$ is singular.] Thus, we may relax the assumption about nonsingularity of $V$ and arrive at the following result (cf. remark 5.4).{\bf Theorem 6.1.} {\it The estimates of theorem 5.1 apply to the parallel arithmetic complexity of computing $\det V$ (with double number of random parameters and double failure probability) and $V^{-1}$ for an $n \times n$ Toeplitz-like matrix $V$. Over the fields of characteristics $p=0$ and $p>n$, the estimates are deterministic. Over the fields of positive characteristics $p \leq n$, the estimates are randomized; they are of Monte Carlo type for testing if $V$ is singular, and they are of Las Vegas type for testing if $V$ is nonsingular. If $V$ is nonsingular,   they are of Las Vegas type also for computing $V^{-1}$ and $\det V$.} {\bf Remark 6.1.} If the inverse of a nonsingular Toeplitz-like matrix $V$ is sought in the ring of integers modulo $p^s$ for a prime $p$, then one may first compute $V^{-1} \bmod p$, performing all the field operations modulo $p$ (with a lower precision), and then lift the solution to $V^{-1} \bmod p^s$ by means of the algorithm of [MC79], which only involves matrix multiplications and subtractions. Moreover, as this was pointed out in [P92a] (cf. also algorithm 2.1 of [P90/92] for $\lambda=p$), all the auxiliary matrices can be computed in the algorithm of [MC79] with their short displacement generators, which enables us to decrease the computational cost dramatically.  {\bf Remark 6.2.} The algorithms of the next section give us a Las Vegas randomized singularity test for a Toeplitz-like matrix. An alternative approach is via recursive triangular factorization of $V$, which also gives us randomized algorithms for computing matrix ranks and determinants (cf. [P85], [P87], [BP94], pages 99, 100, 212, 323, and [P,a]). \noindent {\bf 7. Computation of the rank and the characteristic polynomial of a Toeplitz-like matrix.} In this section we will extend algorithm 5.1 to the computation of the rank of an $n \times n$  Toeplitz-like  matrix $V$ over  any field ${\bf F}$ of any characteristic $p$ and the characteristic polynomial of $V$ provided that $p=0$ or $p>n$. We will write $\hbox{RANK}(n)$ and $\hbox{CHAR}\cdot \hbox{POL}(n)$ to denote these 2 computational problems. For $p=0$ and for $p>n$, the main algorithm of [P90/92]  computes, at the deterministic cost bounded by $C_n$ of (1.2), the characteristic polynomials $C_V$ of $V$ (thus solving $\hbox{CHAR}\cdot\hbox{POL}(n)$  for $p=0$ and $p>n$) and $C_{\tilde{U}}(x)$ of the $2n \times 2n$ symmetric Toeplitz-like matrix  $\tilde{U}=\left(\begin{array}{lr} O & V\\V^T & O\end{array} \right)$. Recall that, for $p=0$, rank $\tilde{U} =2$ rank $V=2n-q$, where $x^q$ is the highest power of $x$ that divides $C_{\tilde{U}}(x)$, and output rank $V=n-(q/2)$, thus solving $\hbox{RANK}(n)$ for $p=0$. It remains to solve RANK($n$) for $p>0$. Our algorithm 7.1 will give us a solution for $p>n$ via a randomized reduction to  $\hbox{CHAR}\cdot \hbox{POL}(n)$. To solve $\hbox{RANK}(n)$ for $0 < p \leq n$, we will shift from the input matrix $V$ to another $n \times n$ Toeplitz-like matrix $W$ having the same rank and having a characteristic polynomial whose associated ratio $B_1(y)/B_0(y)$ [cf. (4.1)--(4.2)] satisfies certain degeneracy restriction (to be specified in definition 7.2). The problem of computing the characteristic polynomials of matrices $W$ of such a class will be referred to as $\hbox{CHAR} \cdot \hbox{POL}^*(n)$. We will show reductions of $\hbox{RANK}(n)$ to $\hbox{CHAR} \cdot \hbox{POL}^*(n)$ (in algorithm 7.1)  and of $\hbox{CHAR} \cdot \hbox{POL}^*(n)$ to $\hbox{RANK}(d)$, $d=\lfloor n/p \rfloor$ (in algorithm 7.2). Based on these reductions, we will define a {\it bipartite recursive process}, which will enable us to reduce the problems to $\hbox{CHAR} \cdot \hbox{POL}^*(k)$ for $k<p$ in $\lceil (\log n)/\log p\rceil -1$ steps, and then the main algorithm of [P90/92] will give us the desired solution. We will start with the following definition and auxiliary results. {\bf Definition 7.1.} A matrix $A$ has a rank $r$ if the maximum size of its nonsingular submatrix is $r \times r$. {\bf Fact 7.1.} {\it Let $A$ be an $n \times n$ matrix having a rank $r$ and let $L$ and $U^T$ be 2 generic $n \times n$ unit lower triangular Toeplitz matrices. Then the $i \times i$ leading principal submatrices of the matrix $UAL$ and the $i \times i$ trailing principal submatrices of the matrix $LAU$ are nonsingular for $i=1, 2, \ldots, r$.} {\bf Proof} was given in [KS91] (and is also available on page 206 of [BP94]) for the first statement of fact 7.1, claiming nonsingularity of the leading principal submatrices of $UAL$. Essentially the same proof also applies to the second statement of fact 7.1, claiming nonsingularity of the trailing principal submatrices of $LAU$. Alternatively, we may apply the first statement of fact 7.1 to $A, \  L$ and $U$ replaced by $\tilde{A}, \  \tilde{L}$ and $\tilde{U}$, respectively; write $L=J\tilde{U}J, \  A=J\tilde{A}J, \  U=J\tilde{L}J$, for $J$ of definition 2.1, $J^2=I$, and then immediately deduce the second statement of fact 7.1.   \hfill  $\Box${\bf Fact 7.2.} {\it Let, under the assumptions of fact 7.1, the $2n-2$ indeterminates-entries of the matrices $L$ and $U$ be replaced by $2n-2$ random parameters chosen from a finite set $S$ of a cardinality $|S|$ (independently of each other and under the uniform probability distribution on $S$). Then, a) the $r \times r$ trailing principal submatrix of $LAU$ is nonsingular with a probability at least $1-2r/|S|$, whereas b) all the $i \times i$ trailing principal submatrices of $LAU$ for $i\leq r$ are nonsingular with a probability at least $1- (r+1)r/|S|$.} {\bf Proof.} Observe that the entries of $LAU$ are quadratic polynomials in the entries of $L$ and $U$. Therefore, the determinant of the trailing principal submatrix of $LAU$ is a polynomial of a degree at most $2i$ in the entries of $L$ and $U$. Combine this observation with facts 2.8 and 7.1 and obtain fact 7.2. [To deduce part b), observe that  $\prod_{i=1}^n (1-2i/|S|) \geq 1-2\sum_{i=1}^r i/|S| =1-(r+1)r/|S|$.]  \hfill $\Box$By definition 7.1, all the $i \times i$ submatrices of $A$ are singular for $i>r$, and due to this property and fact 7.2, we may devise a randomized algorithm for computing $r = {\rm rank} \  A$ if we apply binary search for the maximum $i$ for which the $i \times i$ trailing principal submatrix of $LAU$ is nonsingular. Over any field, we may test nonsingularity of such submatrices of a $k \times k$ Toeplitz-like matrix $A$ by applying our randomized algorithm supporting theorem 6.1.This will give us rank $A$ in $\lceil \log n \rceil $ such successive nonsingularity tests, at the price of increasing the time bounds of theorems 5.1 and 6.1 and decreasing their processor bounds by the same factor  $\lceil \log n \rceil $. In this section, we will achieve the same goal by using Las Vegas randomized algorithms, with no increase of the time bound for $p>n$ and with its increase by factor $\gamma(n,p)$ of (5.5) for $2\leq p \leq n$. As a by-product, we will also solve $\hbox{CHAR}\cdot \hbox{POL}^*(n)$.%{\bf Proof.}  The next randomized algorithm reduces $\hbox{RANK}(n)$ to $\hbox{CHAR}\cdot\hbox{POL}(n)$, over any field ${\bf F}$ of any positive characteristic $p$. {\bf Algorithm 7.1} [extension from $\hbox{CHAR}\cdot\hbox{POL}(n)$ to $\hbox{RANK}(n)$]. {\bf Input:} a field ${\bf F}$ having a positive characteristic $p$ and a cardinality $|{\bf F}|$, a black box subroutine for $\hbox{CHAR}\cdot\hbox{POL}(n)$,  a displacement generator of a length $\ell$ for an $n \times n$ Toeplitz-like matrix $V$, and a finite set $S$ of cardinality $|S|$, consisting of the elements of the field ${\bf F}$ or of its fixed algebraic extension field ${\bf E}$. {\bf Output:} either FAILURE (with a probability bounded by  some fixed upper bound) or $r= {\rm rank} \  V$. {\bf Computations.}  \begin{enumerate} \itemChoose $(8\ell +66)n-2$ random values from the set $S$; $2n -2$ of them are allocated for the entries of 2 random unit $n \times n$ lower triangular Toeplitz matrices $L$ and $U^T$, $2(2\ell +14)n$ of them form a pair of $n \times (2\ell +14)$ random matrices $X$ and $Y^T$, and $2(2\ell +18)n$ form a pair of $n \times (2\ell +18)$ random matrices $X_+$ and $Y_+^T$.  \item %No 2 Compute displacement generators for 2 matrices: of a length at most  $\ell +4$ for $LVU$ and of a length at most $5\ell +39$ for the matrix $$ W=GLVUH \ ,  \eqno (7.1) $$where $G-ZGZ^T = XY^T$ and $H-ZHZ^T = X_+Y_+^T$ (compare definition 2.2).\item %no 3 Apply the given subroutine for $\hbox{CHAR}\cdot\hbox{POL}(n)$ in order to compute the characteristic polynomial $C_W(x)$ of the matrix $W$; then compute the maximum $q$ such that $x^q$ divides $C_W(x)$ and then the value $r_-=n-q$. \item %no 4 Test whether the $r_-\times r_-$ trailing principal submatrix of $LVU$ is nonsingular, by applying the algorithm of section 6 that supports theorem 6.1. If so, output  $r=r_-$; otherwise, output FAILURE. \end{enumerate}{\bf Remark 7.1.} Unlike the singularity test of theorem 6.1, we have Las Vegas randomization in algorithm 7.1: unless FAILURE is output (with a low probability), the $r_- \times r_-$ submatrix of $LVU$ is nonsingular, and this enables us to verify that $r_-=r$. We will first show correctness of algorithm 7.1 assuming that the unit lower triangular Toeplitz matrices $L$ and $U^T$, as well as the matrices $X$, $Y$, $X_+$, and $Y_+$, are generic matrices, filled with $(8\ell +66)n-2$ indeterminates, rather than with $(8\ell +66)n-2$ random elements from $S$. In this case, clearly, we have $ r = {\rm rank} \   V = {\rm rank} \ W \ , % \eqno (7.6) $ for $W$ of (7.1), and furthermore, $r_- \leq r = {\rm rank} \ W$. Therefore, it suffices to show an assignment for the entries of $X$, $Y$, $X_+$, and $Y_+$ for which $r_-=r  = {\rm rank} \  W$, and then we may conclude that $r=r_-$ for generic $X$, $Y$, $X_+$, and $Y_+$.  We first write $$ LVU = \left( \begin{array} {cc}  B & C  \\ D &  E  \end{array} \right), $$ where $E$ is an $r\times r$ matrix. Due to fact 7.1, $E$ is nonsingular. Now, we write $ H=\tilde{H}\tilde{F}$,  $$G = \left( \begin{array}{cc} I& -CE^{-1} \\O&I\end{array} \right) \ , \ \  \tilde{H}= \left( \begin{array}{cc} I& O \\-E^{-1}D & E^{-1}\end{array} \right) \ , \ \ \tilde{F}= \left( \begin{array}{cc} O& O \\O& F_r\end{array} \right) \ , %\eqno (7.7) $$   where $F_r$ is a generic $r \times r$ Frobenius matrix (compare definition 2.6). Then, since rank $(LVU)=r={\rm rank}\  E$, we have  $GLVU\tilde{H} = \left( \begin{array}{cc} O& O \\O&I_r\end{array} \right)$, $I_r$ being the $r\times r$ identity matrix. Therefore, $W=\tilde{F}$, for $W$ of (7.1), and consequently (cf. definition 7.1), $$C_W(x) = x^{n-r} \sum_{i=0}^r c_{r-i}(W) x^i \ . \eqno (7.2) $$Here, $c_0(W)=1$, and $c_1(W), \ldots, c_r(W)$ are indeterminates. Therefore, $r_-=r$, as was required. Let us next show that also the displacement ranks of the matrices $G$, $H$, and $W$, under (7.1), are bounded as this was required in algorithm 7.1. We will use the following simple result from [P93a] (compare lemma 4.6.2 on page 344 of [BP94]).{\bf Fact 7.3.} {\it Let $$W = \left( \begin{array}{cc} W_{0,0} & W_{0,1}  \\W_{1,0} & W_{1,1} \end{array} \right) \ , \  \  W(0) = \left( \begin{array}{cc} W_{0,0} & O  \\W_{1,0} & O \end{array} \right) \ , \ \  W(1) = \left( \begin{array}{cc} O & W_{0,1}  \\O & W_{1,1} \end{array} \right)  $$and let $W(i,j)$ be obtained from $W$ by substituting $W_{g,\  h}=O$ unless $g=i \ , \  h=j$ for $i ,j = 0, 1 $, so that $W(0,0)={\rm diag}(W_{0,0},O)$, $W(0,1)={\rm diag}(W_{0,1},O)J$, $W(1,0)={\rm diag}(O, W_{1,0})J$, $W(1,1)={\rm diag}(O, W_{1,1})$, $W(j)=W(0,j) +W(1,j), \  j=0, 1$. Let $\rho$, $\rho(j)$, $\rho(i,j)$, and $\rho_{i,j}$ denote the displacement ranks of the matrices $W$, $W(j)$, $W(i,j)$, and $W_{i,j}$, respectively, for $i,j=0, 1$. Then $\rho_{i,j} \leq \rho + i + j$, $\rho(i,j) \leq \rho +2$, $\rho(j) \leq \rho +j$, $\rho(i,j) \leq \rho_{i,j} +2 - i -j$, for $i, j = 0, 1$.}Next, represent $H$ as $H_0H_1\tilde{F}$, where $$H_0 =  \left( \begin{array}{cc} I & O  \\O & E^{-1} \end{array} \right) \ , \ \ H_1 =  \left( \begin{array}{cc} I & O  \\-D & I \end{array} \right) \ . $$ Now, recall definition 2.2, successively apply facts 2.4 and 7.3, and obtain that $\rho_E \leq \rho_{LVU} \leq \rho +4$, $\rho_C \leq \rho_{LVU} +1 \leq \rho +5$, $\rho_D \leq \rho_{LVU} +1 \leq \rho +5$,  $\rho_{H_1} \leq \rho +7$,  $\rho_{E^{-1}} \leq \rho +6$,  $\rho_{-CE^{-1}} \leq 2\rho +12$, $\rho_{H_0} \leq \rho +8$,$\rho_{\tilde{F}} = \rho_{F_r}\leq 2$,  $\rho_G \leq 2\rho +14$,  %%$\ell_{E^{-1}} \leq \ell +6$,  $\rho_H \leq 2\rho + 19$,  $\rho_W \leq \rho_G +\rho_{LVU} +\rho_H +2 \leq 5\rho + 39$.  Note that these bounds on $\rho_G$, $\rho_H$, and $\rho_W$ are consistent with the assumptions made at stages 1 and 2 of algorithm 7.1. This completes the correctness proof for algorithm 7.1 in the case where all the random entries of the matrices $L$, $U$, $X$, $Y$, $X_+$ and $Y_+$ are replaced by indeterminates.Now, let us revisit algorithm 7.1 assuming that the $2n-2$ subdiagonal entries of the unit lower triangular Toeplitz matrices $L$ and $U^T$ and all the entries of the matrices $X$, $Y$, $X_+$, and $Y_+$ are random values (from the set $S$), rather than indeterminates. Such an assumption cannot increase the ranks and the displacement ranks, so we still have $\rho_G \leq 2\rho +14$, $\rho_H \leq 2\rho +19$, $\rho_W \leq 5\rho + 39$, as was required, and also we have rank $W\leq r = {\rm rank}\  V$. On the other hand, the entries of the matrix $W$ are polynomials of degrees at most 6 in the entries of the matrices $L$, $U$, $X$, $Y$, $X_+$, and $Y_+$. Furthermore, $c_r(W)$ is a polynomial of a degree at most $r$ in the entries of $W$. Therefore, due to facts 2.8, 7.2, and definition 7.1, we have $${\rm Probability} \ ({\rm rank} \  W={\rm rank} \ E=r) \geq 1-6r/|S| \ , $$ $${\rm Probability} \ (c_r(W) \neq 0)= \ {\rm Probability} \ (r_- \geq {\rm rank}\  W) \geq 1-6r/|S| \ . $$Since, clearly, $r_- \leq r$, we deduce from the 2 latter estimates that $r_-=r$ with a probability at least $1-\hat{P}_0=(1-6r/|S|)^2 \geq 1-12r/|S|$, that is, we have $$ \hat{P}_0 \leq \frac{12 n}{|S|} \eqno (7.3) $$ since $r\leq n$. By summarizing our analysis and including simple estimates for the {\it computational cost} of performing stages 2, 3, and 4 of algorithm 7.1, based on application of fact 2.4 at stage 2 and of the main algorithm of [P90/92] at stages 2 and 3, we obtain the following theorem.{\bf Theorem 7.1.} {\it Let $\beta$ be a positive scalar, let $V$ be an $n \times n$ Toeplitz-like matrix given with its displacement generator of a length $\ell$, let ${\bf F}$ be a field of a characteristic $p$, where $p=0$ or $p>n$, and let$ N_0 = 2^{\beta}(12n)$. Then one may compute over ${\bf F}$: \\ a) rank $V$ if $p=0$ and the characteristic polynomial $C_V(x)$ if $p=0$ and if $p>n$, both at a deterministic cost bounded by $C_n $ of (1.2), \\ b) rank $V$, at a randomized Las Vegas cost [not including the cost of generating random parameters and of computing the auxiliary algebraic extension field (cf. remark 5.1) and] bounded either by $C_n$ of (1.2) if $p>n$ and if ${\bf F}$ has at least $N_0$ distinct elements (that is, if $|{\bf F}| \geq N_0$) or [compare (5.4)] by $$ C_{n, |{\bf F}|} = O((\log n)^2, (n^2\log \log n)(\log N_0)/((\log n)\log |{\bf F}|)) %\eqno (7.2)$$if $p> n$ and if $|{\bf F}| < N_0$; in the latter case, where $|{\bf F}| < N_0$,  we have  $$C_{n, |{\bf F}|} = O((\log n)^2 ,(n^2 \log\log n)/\log |{\bf F}| )  \    {\rm if} \   \beta=O(\log n)\ .$$ If $p>n$, then,   under both assumptions on $|{\bf F}|$, the failure probability is at most $2^{-\beta}$ and at most $(8\ell + 66)n -2$ random parameters are required.}It remains to consider the case of the fields ${\bf F}$ of positive characteristics $p\leq n$. We will start with defining   $\hbox{CHAR}\cdot\hbox{POL}^*(n)$ and with a reduction of $\hbox{CHAR}\cdot\hbox{POL}^*(n)$ to $\hbox{RANK}(d)$, for $d=\lfloor n/p \rfloor$.Let us represent the characteristic polynomial $C_W(x)$ of the matrix $W$ according to (7.2), where $c_0(W)=1$ and where $c_{r-i}(W)$, for $i=0, 1, \ldots, r-1$, are polynomials of degrees at most $r-i$ in the entries of $W$. Furthermore, let us write $C(x)=C_W(x)$ and $c_i=c_i(W)\ , \   i=0, 1, \ldots, n$, in order to simplify the notation. Once again, we will first assume that $c_1, \ldots, c_r$ are indeterminates (and we will apply fact 2.8 when we come back to randomized $c_1, \ldots, c_r$ later on).  [More precisely, we will temporarily assume that the matrix $W$ satisfying (7.2) is given with its displacement generator filled with indeterminates, and under this assumption, we   will deduce all the desired properties of the coefficients of $C(x)$.] Then, (7.2) implies that  $c_{k+ip}=0$ if and only if $k+ip >r$ or,  equivalently, $$c_{k+ip} =0 \ {\rm if \  and \  only \   if } \  i>\lfloor (r-k)/p \rfloor \ , \eqno (7.4) $$for $k=0, 1, \ldots, p-1$. We recall the notation of section 4 [see (4.1)--(4.3)] and deduce from (7.4) that $$ B_k(y) = 0 \   {\rm (identically \  in } \  y) \  {\rm if \   and \  only \  if} \   k>r \ , \eqno (7.5) $$ $$ \deg B_k(y) = \lfloor (r-k)/p\rfloor \  {\rm otherwise}\ ,\eqno (7.6) $$for $k=0, 1, \ldots, p-1$. {\bf Definition 7.2.}  $\hbox{CHAR}\cdot\hbox{POL}^*(n)$ is the problem of computing the characteristic polynomial of an $n \times n$ Toeplitz-like matrix $W$ (having an unknown rank $r$) for which simultaneously (7.2) holds, the associated polynomials $B_1(y)$ and $B_0(y)$ have only constant common divisors, and $B_i(y)$ satisfy (7.5) and (7.6), for $i=0$ and/or $i=1$. Clearly, the assumptions of definition 7.2 hold if $c_1, \ldots, c_r$ are indeterminates, under (7.2). Due to corollary 2.2, we also have {\bf Fact 7.4.} {\it Under (7.2), the assumptions of definition 7.2 holds if and only if the associated polynomials $B_0(y)$ and $B_1(y)$ define a nonsingular $s \times s$ Sylvester (resultant matrix), for $s= \lfloor \frac{r}{p}\rfloor + \lfloor \frac{r-1}{p}\rfloor$.} Now, we are ready to present the desired algorithm (which later on we will combine with algorithm 7.1 for $\tilde {\ell} = 5\ell +39$). {\bf Algorithm 7.2} [extension from $\hbox{RANK}(d)$ to $\hbox{CHAR}\cdot\hbox{POL}^*(n)$ ]. {\bf Input:} a field ${\bf F}$ of a positive characteristic $p\leq n$;  a set $S$ (as in algorithm 7.1), a displacement generator of a length $\tilde{\ell} \geq 2$ for an $n \times n$ Toeplitz-like matrix $W$ satisfying the assumptions of definition 7.2, and a black box subroutine for $\hbox{RANK}(d)$, $d=\lfloor n/p\rfloor$.   {\bf Output:}  $C_W(x)$.{\bf Computations.} \begin{enumerate}\item %no 1 Apply the main algorithm of [P90/92] to compute [at deterministic cost $C_n$ of (1.2)] the values $s_k=s_k(W)={\rm trace} \  (W^k),\  k=1, 2, \ldots, 2n$. (Note that $s_k$ is the $k$th power sum of the zeros of $C(x)=C_W(x)$ for $k=1, 2,\ldots, 2n$.) \item %no 2 Compute (at a deterministic cost bounded according to theorem 4.1) the coefficients $a_1, a_2, \ldots, a_m$, $m=n+\lfloor (n-1)/(p-1)\rfloor < 2n$, of the formal power series $A(x)$ satisfying (4.5) and (4.7). \item %no 3 Compute the natural $d_k$, $k=0, 1$, from (4.1)--(4.3); then apply fact 2.10, the given subroutine for $\hbox{RANK}(d)$,  and algorithm 5.1 for $n$ replaced by $d = \lfloor n/p \rfloor$ in order to compute the $(d_1, d)$ entry $(B_1(y),\  B_0(y))$ of the Pad\'{e} approximation table for the formal power series $A_1(y)$ defined by $A(x)$ [see (4.4), (4.6)].\item %no 5Compute the coefficients of the polynomials $B_k(y)=B_0(y)A_k(y) \bmod y^{d_k+1}$, for $k=2, 3, \ldots, p-1$, which gives us all the coefficients of $B(x)$ and $C(x)$.  Output $C_W(x)=C(x)$. \end{enumerate}Under the assumption that (7.2) holds for the input matrix $W$ defining $\hbox{CHAR} \cdot \hbox{POL}^*(n)$, {\it correctness} of algorithm 7.2 immediately follows from (7.4)--(7.6) and theorems 4.1 and 5.1. Now, we shall come back to the assumption that the given displacement generator for the matrix $W$ is randomized (rather than generic), that is, filled with random elements of a set $S$ fixed in the fields ${\bf F}$ or ${\bf E}$ (rather than with indeterminates). Then, the assumptions of definition 7.2 will hold with a high probability, and application of theorems 4.1 and 5.1 will enable us to bound the computational cost of performing algorithm 7.2.  To solve both problems $\hbox{RANK}(n)$ and $\hbox{CHAR}\cdot\hbox{POL}^*(n)$, where $0< p \leq n$, we will recursively apply algorithms 7.1 and 7.2 (by using them as black box subroutines for each other) until the initial input size $n$ for $\hbox{CHAR}\cdot\hbox{POL}^*(n)$ decreases below $p$, which will occur in less than $(\log n)/\log p$ bipartite recursive steps. Then a single application of the main algorithm of [P90/92] will complete the solution. By summarizing the estimates for the computational cost, the number of random parameters involved, and the failure probability, we will arrive at the next theorem (cf. remark 5.4). {\bf Theorem 7.2.} {\it Let $\beta$ be a positive scalar, let $V$ be an $n \times n$ Toeplitz-like matrix given with its displacement generator of a length $\ell$, let ${\bf F}$ be a field of a characteristic $p$, $0<p\leq n$, and let $$ N_1=\lceil 2^{\beta +1}n\left( \frac{2(p^3 +1)n^2}{(p^3 -1)^2} + \frac{6n}{p^2-1} + \frac{7p-5}{(p-1)^2} p \right) \rceil \ .$$Then one may compute rank $V$ at a randomized Las Vegas cost [not including the cost of generating random values and of computing the auxiliary algebraic extension field and]  bounded either by $$ C_{n,p}^* = O((\gamma(n,p) \log n )^2, ((n/\gamma(n,p))^2(\log \log n)/\log n))\eqno (7.7) $$ if $|{\bf F}| \geq N_1$ or by $$ C_{n, p, |{\bf F}|}^* =  O((\gamma(n,p) \log n )^2, ((n/\gamma(n,p))^2(\log N_1) (\log \log n)/((\log n)(\log {\bf F})))\eqno (7.8) $$ if $|{\bf F}| < N_1$. In both cases, $\gamma(n,p)$ is defined by (5.5), a failure probability is  bounded from above by $2^{-\beta}$, and at most $(10\ell +98+ \frac{146}{p-1})n -2\lfloor (\log n)/\log p\rfloor \leq (10\ell+244)n$ random parameters are required, for $p\geq 2$, $n\geq 1$.} {\bf Proof.} The cost of our computation of rank $V$ is clearly dominated by the overall cost of performing algorithm 5.1 at stage 4 of algorithm 7.1 and at stage 3 of algorithm 7.2, in all the recursive steps, which we define as follows. At the initial, 0-th step, algorithm 5.1 is applied to the $n \times n$ Toeplitz-like matrix $LVU$ [see (7.1) and stage 4 of algorithm 7.1] , whose displacement generator of a length $\ell_0 \leq \ell +4$ ia available . At the $i$-th recursive step, $i=1, 2, \ldots, h$, we apply algorithm 5.1 at stage 4 of algorithm 7.1 to an $n_i \times n_i$ matrix defined with its displacement generator of a length $\tilde{\ell}_i \leq 6$, similarly to the $LVU$, and at stage 3 of algorithm 7.2 to an $n_i \times n_i$ Toeplitz matrix (given with their displacement generator of a length $\ell_i \leq 2$), where  $$ n_0=n, \  \  n_i \leq \lfloor n_{i-1}/p \rfloor, \  \   i=1, \ldots, h; \   \   h \leq \lfloor (\log n) /\log p\rfloor -1 \ . \eqno (7.9) $$ Then, application of (5.3), (5.4), and Brent's principle yields (7.7) and (7.8). Note that in addition to the inner loop of bipartite recursion in algorithm 5.1, we have the outer loop of the bipartite recursive transition $\hbox{RANK}(n) \longrightarrow \hbox{CHAR}\cdot \hbox{POL}^*(n) \longrightarrow \hbox{RANK}(d)$, which leads to increasing [by factor $\gamma(n,p)$ of (5.5)] the time bound of theorem 5.1. We will next estimate the number of random parameters involved in the bipartite recursive process of application of algorithms 7.1 and 7.2 and the failure probability.Due to theorem 5.1, recursive application of algorithm 5.1 at stages 3 and 4 of algorithm 7.2 contributes at most $(2\ell_0+24)n + \sum_{i=1}^h(2\ell_i +2\tilde{\ell}_i +48)n_i$ extra random parameters,  where $h$ and $n_0, \ldots, n_h$ are defined by (7.9) and where $\ell_0=\ell +4, \  \ell_i=2, \  \tilde{\ell}_i = 6$ for $i>0$ (see above). The total number of these random parameters is thus bounded by $(2\ell +32)n + \sum_{i=1}^h 64n/p^i < (2\ell +32 +64/(p-1))n$. Besides that, we generate $(8\ell_i + 66)n_i - 2$ random values at $i$-th recursive application of stage 1 of algorithm 7.1, $i=0, 1, \ldots, h$. This gives us $ \sum_{i=0}^h ((8\ell_i +66)n_i -2) = (8\ell +66)n + 82 \sum_{i=1}^h n_i - 2 \lfloor (\log n)/\log p \rfloor < (8\ell +66)n + 82 \frac{n}{p-1} - 2 \lfloor (\log n)/\log p \rfloor $ additional random values. By summing both bounds, we obtain the bound $(10\ell +98 + \frac{146}{p-1})n - 2 \lfloor (\log n)/\log p \rfloor $ of theorem 7.2 on the overall number of random parameters involved.  Let us next estimate the failure probability. Due to theorem 5.1, all the applications of algorithm 5.1 during the recursive process  contributea chance of failure with a probability at most $ P^* = 1 - (1- \frac{2np}{N_1}(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)})) \prod_{i=1}^h (1-\frac{2n_ip}{N_1}(\frac{1}{p-1} + \frac{2n_i^2}{(p^3-1)}))^2 \leq \frac{2np}{N_1}(\frac{1}{p-1} + \frac{2n^2}{(p^3-1)}) + \sum_{i=1}^h \frac{4n_ip}{N_1}(\frac{1}{p-1} + \frac{2n_i^2}{(p^3-1)}), \  {\rm for} \  h, n_1, \ldots, n_h$ of (7.9). It follows that $$ P^* < \frac{2np}{N_1}((\frac{(p+1)}{(p-1)^2}+ 2(p^3 +1) (\frac{n}{p^3-1})^2) \ .  \eqno (7.10) $$Likewise, we extend our earlier bound (7.3) so as to cover the contribution of the $h+1$ recursive calls for algorithm 7.1. This gives us $ \hat{P}=1-\prod_{i=0}^h (1-\hat{P}_i) \leq \sum_{i=0}^h \hat{P}_i$, where $P_i = \frac{12n_i}{N_1}$, and $h$, $n_0, \ldots, n_h$ are defined by (7.9). Therefore, the contribution $\hat{P}$ to the overall probability of failure is bounded by $$ \hat{P} \leq  \frac{12np}{N_1(p-1)} \ . \eqno (7.11) $$Finally, we shall estimate the contribution to the failure due to the possible degeneration of some entries of the Pad\'{e} tables that we need to compute in the recursive applications of algorithm 7.2. Due to fact 7.4, degeneration of such an entry is equivalent to vanishing of the associated resultant, that is, of the determinant of the associated Sylvester matrix $R$. In the generic case, $R$ is nonsingular. In the randomized case, we will apply fact 2.8 to estimate the failure probability. The resultant, $\det R$, is a polynomial in the coefficients $c_j$ of the characteristic polynomial, and $\det R$ has a degree less than $2n_i+1$ in these coefficients $c_j$ at the $i$-th recursive step, $i=1,2, \ldots, h$.Since $c_j$ is a polynomial of a degree at most $j$ in the entries of $W$, $j=1,2, \ldots, n$, and since the entires of $W$ are polynomials of degrees at most 6 in the entries of $L$, $U$, $X$, $Y$, $X_+$, and $Y_+$, we obtain that, at the $i$-th recursive application  of algorithms 7.2, the Pad\'{e} entry degenerates with a probability at most $12n^2/N_1$, where $N_1=|S|$. Therefore, we will avoid degeneracy, at all recursive applications of algorithm 7.2,with a probability at least $$ 1-\tilde{P} = \prod_{i=1}^n (1-\frac{12n_i^2}{N_1})\geq 1-\frac{12}{N_1} \sum_{i=1}^h n_i^2 \geq1-\frac{12n^2}{N_1}\sum_{i=1}^h n_i^2 $$ $$ \  \  \geq 1-\frac{12n^2}{N_1}\sum_{i=1}^h p^{-2i} \geq  1- \frac{12n^2}{(p^2-1)N_1} \ .$$ Combining the latter bound with (7.10) and (7.11) gives us the upper bound $$ P^* + \hat{P} + \tilde{P} < \frac{2n}{N_1}(\frac{(p+1)p}{(p-1)^2}+ 2(p^3+1) (\frac{n}{p^3-1})^2p + \frac{6p}{p-1}+  \frac{6n}{p^2-1}) $$ $$ \  \  \   =\frac{2n}{N_1}(2(p^3 +1) (\frac{n}{p^3-1})^2p+ \frac{6n}{p^2-1}+\frac{7p -5}{(p-1)^2} p)$$  %%$$ \  \  \   \   \   \  \leq \frac{2n}{N_1}(\frac{9n^2}{7} %%+ \frac{3}{n}\lfloor \log_2 n \rfloor + 13 + \frac{6n}{N_1}) %\ . \eqno (7.18) %%$$on the overall probability of failure in the recursive application of algorithms 7.1 and 7.2.Therefore, in order to keep the overall probability of failure of algorithm 7.1 below a fixed positive $\beta$, it suffices to choose $N_1 = |S|$ satisfying $N_1 \geq  2^{\beta+1}n(2(p^3+1)(\frac{n}{p^3-1})^2 + \frac{6n}{p^2-1}+\frac{7p -5}{(p-1)^2} p)$. %%\frac{(p+1)p}{(p-1)^2} + %%\frac{4n^3(p^3 +1)}{(p^3-1)^2 p^2} %%+ 6 \lfloor \frac{\log n}{\log p} \rfloor + %%18\frac{np}{p-1}  \}$. This completes the proof of theorem 7.2. \hfill $\Box${\bf Remark 7.2.} The results of theorems 7.1 and 7.2 can be immediately extended to any $k \times h$ Toeplitz-like matrix $T$, where $k+h \leq n$, since we may embed the matrix $T$ into an $n \times n$ Toeplitz-like matrix $V$ obtained by appending zero rows and/or columns to $T$.{\bf Remark 7.3.} To extend the proposed algorithm to the computation (over any field) of the characteristic polynomial of a fixed input matrix $V$, we may follow the line of remark 5.3 and replace $W$ and all other auxiliary random matrices by the identity matrices. Then the algorithms of this section will solve $\hbox{I}\cdot \hbox{P}\cdot \hbox{SUMS}(n)$ and $\hbox{CHAR}\cdot \hbox{POL}(n)$ unless some auxiliary Pad\'{e} entries degenerate. The correctness test (verifying that we have arrived at $C_V(x)$ indeed) is simple, since it reduces to solving a triangular Toeplitz system of Newton's identities [see (2.4), (2.5)]. For more details and for application to the algebraic coding theory, we refer the reader to [P,b]. \vspace{.2in} {\bf 8. The null-space of a Toeplitz-like matrix.}  {\bf Definition 8.1} (compare [GL89] or [BP94]). The null-space ${\bf N}$(A) of a matrix $A$ is the linear space of all vectors $\vec{x}$ satisfying $A\vec{x} = \vec{0}$ [so that the 2 vector equations $A\vec{x} = \vec{f}$ and $A\vec{y} = \vec{f}$together imply that $\vec{x} - \vec{y} \in {\bf N}$(A)]. Next, we will consider computation of the null-space ${\bf N}$(V) for a Toeplitz-like matrix $V$, over any field ${\bf F}$. We will apply the following algorithm (compare [KP91]). {\bf Algorithm 8.1.}  {\bf Input:} a field ${\bf F}$, a finite set $S$ in ${\bf F}$ or in its algebraic extension field, and a displacement generator of a length $\ell= O(1)$ for an $n \times n$ Toeplitz-like matrix $V$. {\bf Output:} a basis for the null-space ${\bf N}$(V) over ${\bf F}$. {\bf Computations.} \\ 1. Apply the bipartite recursive algorithm of section 7 in order to compute $r={\rm rank}\ V$ and a displacement generator of a length at most $\ell +4$ for the matrix $LVU$, where $L$ and $U$ are 2 random matrices defined in fact 7.2. In particular, at stage 4 of the first application of algorithm 7.1 in the outer loop of this bipartite recursive algorithm, the $r \times r$ trailing principal submatrix $E$ of the matrix  $LVU$ of (7.1) is computed, and the algorithm supporting theorem 6.1 is applied in order to test if $E$ is nonsingular. If so, the algorithm computes a displacement generator of a length at most $\ell +6$ for $E^{-1}$; if $E$ is singular, the algorithm outputs FAILURE. \\2. Write  $LVU=\left( \begin{array}{lr} B & C\\D & E \end{array} \right)$ and  $G=\left( \begin{array}{lr} I & O\\-E^{-1}D & I \end{array} \right)$;  compute a displacement generator of a length at most $2\ell +4$ for the matrix $-E^{-1}D$ and a displacement generator of a length at most $2\ell +6$ for the matrix $G$  (compare facts 2.4 and 7.3). \\ 3. Compute the matrix $M=UG\left( \begin{array}{c} I \\O \end{array} \right)$, where $I$ is $(n-r) \times (n-r)$ identity matrix. If $VM=O$, output the $n-r$ columns of $M$ as a basis for the null-space of $V$; otherwise, output FAILURE. To verify {\it correctness} of this randomized Las Vegas algorithm, observe that $LVUG=\left( \begin{array}{lr} O & C\\O & E \end{array} \right)$ for rank $(LVUG) \leq r = \ {\rm rank} \  E$.  Due to fact 2.4, the computational cost of performing algorithm 8.1 is dominated by the cost of performing its stages 1 and 2. The latter cost is bounded according to theorems 6.1, 7.1, and 7.2 and fact 2.4. The computational cost of solving $\hbox{RANK}(n)$ at stage 1 (see theorems 7.1 and 7.2) dominates the overall cost of computations by algorithm 8.1, though $2n-2$ random values are required at stage 1 for generating the matrices $L$ and $U$ even where $p=0$, and this may introduce  some additional chance of failure, with a probability at most $2r/|S|$. Furthermore, over the fields of positive characteristics $p\leq r$, the inversion of the $r \times r$ Toeplitz-like matrix $E$, whoise displacement generator of a length at most $\ell +6$ is readily available (cf. facts 2.4 and 7.3), requires $(2(\ell+6) +24)r= (2\ell+36) r \leq (2\ell+36)n$ additional random parameters and contributes additional chance of failure, with a  probability at most $\frac{2rp}{|S|} (\frac{1}{p-1} + \frac{2r^2}{p^3 -1})$.By summarizing all these observations, we obtain the following result (cf. remark 5.4). {\bf Theorem 8.1.} {\it The estimates for the computational cost, for the number or random parameters and for the failure probability of theorems 7.1 and 7.2  can be extended to the problem of computation of the null-space of $n \times n$ Toeplitz-like matrix having an unknown rank $r$ and given with its displacement generator of a length $\ell$. In this extension, $(2\ell +36)r$ additional random parameters are required, the failure probability increases by at most $\frac{2r}{|S|} ( \frac{1}{p-1} + \frac{2r^2}{p^3-1})$, over the field of positive characteristics $p\leq r$, and furthermore,  the solution algorithm is Las Vegas randomized  even for $p=0$, in which case the algorithm uses $2n -2$ random parameters and has a failure probability at most $2r/|S|$.}\vspace{.2in} {\bf 9. Solution of singular Toeplitz-like linear systems over an arbitrary field.} For a singular linear system $$ V\vec{x}=\vec{f} \ , \eqno (9.1) $$ we shall determine if the system is consistent and, if so, shall seek its solution vectors $\vec{x}$, which are all represented as $$ \vec{x}=\vec{y} + \vec{z} \ , \eqno (9.2) $$ where $\vec{y}$ is any fixed specific solution and $\vec{z}$ is any vector from ${\bf N}$(V), the null-space of $V$ (cf. definition 8.1). The algorithm of the previous section can be applied in order to compute a basis of vectors for ${\bf N}$(V). In this section, we will seek a specific solution to the system (9.1) or indication that this system is inconsistent.The linear system $V\vec{x}=\vec{f}$ of (9.1), where $V$ is a $k \times h$ Teoplitz-like matrix, can be immediately reduced to the linear system  $\left( \begin{array}{lr} O & V \\ V^T & O \end{array} \right) \left(\begin{array}{c}\vec{z}\\ \vec{x} \end{array}\right) = \left(\begin{array}{c}\vec{0}\\ \vec{f} \end{array}\right) $, %%$W\vec{y}=(\vec{0}^T, \vec{f}^T)^T$, where $W$ is the with an $n \times n$ Toeplitz-like coefficient matrix for  $n=k+h$. We will assume that already the input system (9.1) has $n$ equations and $n$ unknowns, and we will apply the following randomized Las Vegas algorithm. {\bf Algorithm 9.1.} {\bf Input:} a field ${\bf F}$ having a characteristic $p$, a displacement generator of a length $\ell = O(1)$ for an $n \times n$ Toeplitz-like matrix $V$, an $n$-dimensional vector $\vec{f}$, and a finite set $S$ of elements of ${\bf F}$ or of its algebraic extension field. {\bf Output:} either FAILURE, or a vector $\vec{x}$ satisfying the vector equation $V\vec{x}=\vec{f}$ of (9.1), or INCONSISTENT, indicating inconsistency of the system (9.1). {\bf Computations.}  \begin{enumerate} \item %no 1 Stage 1 is as in algorithm 8.1.\item %no 2 Note that the system $V\vec{x}=\vec{f}$ of (9.1) is equivalent to the linear system $$ LVU\vec{z}=L\vec{f} \ , \eqno (9.3) $$ provided that $\vec{x}=U\vec{z}$. We will seek a specific solution to (9.3) in the form $\vec{z}=\left( \begin{array}{c} \vec{0}\\ \vec{y}\end{array} \right)$, where $\vec{y}$ is an $r$-dimensional vector. Compute the vector $L\vec{f}$, let $\vec{h}$ denote the $r$-dimensional trailing subvector of the vector $L\vec{f}$. \item %no 3  Multiply $E^{-1}$ by $\vec{h}$ in order to compute the solution vector $\vec{y}=E^{-1}\vec{h}$ to the nonsingular Toeplitz-like linear system of $r$ equations, $E\vec{y}= \vec{h}$ (or to output FAILURE, with a low probability). \item %no 4Compute the vector $\vec{x}=U\vec{z}$, where$\vec{z} = \left( \begin{array}{c} \vec{0}\\ \vec{y}\end{array} \right)$.\item %no 5Compute the vector $V\vec{x}$. If $V\vec{x}= \vec{f}$, output $\vec{x}$. Otherwise, output INCONSISTENT. \end{enumerate}To verify {\it correctness} of algorithm 9.1, rewrite (9.3) as $\left( \begin{array}{lr} B & C\\ D & E \end{array} \right) \  \left( \begin{array}{c} \vec{0}\\ \vec{y} \end{array} \right) \  = \  \left( \begin{array}{c} \vec{g}\\ \vec{h} \end{array} \right)$,   where $\left( \begin{array}{c} \vec{g}\\ \vec{h} \end{array} \right) = \vec{f}$,  and deduce that $E\vec{y} = \vec{h}$. Now, inspection (with application of fact 7.2) immediately shows {\it correctness} of algorithm 9.1, which outputs INCONSISTENT only if the linear system $V\vec{x}= \vec{f}$ is inconsistent indeed.Due to facts 2.3, 2.4, and 7.3, stages 2, 3, 4, and 5 can be performed at the deterministic cost $O(\log n, n\log \log n)$. Stage 1 and the cost of its performance are the same as for algorithm 8.1, and we arrive at the following result.  %$(2 \tilde{\ell}+24)r=(2\ell +32)r$ {\bf Theorem 9.1.} {\it For a Toeplitz-like linear system of $k$ equations and with $h$ unknowns and for $n=k+h$, one may either detect inconsistency of this system or compute its solution (over a fixed field ${\bf F}$ of constants having a characteristic $p$)  by applying a randomized Las Vegas algorithm that runs at the cost bounded according to the estimates of theorem 8.1. The estimates of theorem 8.1 for the number of random parameters used and failure probability are also applied.} \vskip2mm\noindent{\bf 10. Least-square solutions to a singular Toeplitz-like system.}In this section, $\vec{f}, \  \vec{x}, \  A$ and $V$ will denote vectors and matrices filled with real or rational entries, which means, in particular, that we will deal with a field of constants having characteristic $p=0$. One may immediately extend our study to the case of complex vectors and matrices by replacing transposes of matrices by their Hermitian transposes (see the definition of Hermitian transposes in [GL89] or [BP94]). {\bf Definition 10.1} ([GL89], page 243, or [BP94], page 90). For an $k \times h$ matrix $A$, its Moore-Penrose {\it generalized (or pseudo) inverse} $X=A^+$ is a unique $h \times k$ matrix satisfying the four Moore-Penrose conditions: $$AXA=A, \  \   XAX = X, \  \   (AX)^T=AX, \  \   (XA)^T=XA \ .$$Clearly, $A^+=A^{-1}$ for a nonsingular matrix $A$. {\bf Definition 10.2.} For a vector $\vec{v}=[v_0, \ldots, v_{n-1}]^T$, its {\it Euclidean norm} equals $\| \vec{v} \|_2 =(\sum_{i=0}^{n-1} |v_i|^2)^{1/2}$. {\bf Definition 10.3.} A {\it least-squares solution} to a linear system of equations $V\vec{x}=\vec{f}$ of (9.1), where $V$ can be a singular matrix, is given by a vector $\vec{y}$ that minimizes the  Euclidean residual norm $\|V\vec{x} -\vec{f}\|_2$ over all candidate vectors $\vec{x}$, that is, $\|V\vec{y} -\vec{f}\|_2 =\min_{\vec{x}} \|V\vec{x} -\vec{f}\|_2$.  A least-squares solution is generally nonunique [cf. (9.2)]. {\bf Fact 10.1} ([GL89], page 242). {\it $\vec{y}=V^+\vec{f}$ is a least-squares solution to a linear system $V\vec{x}=\vec{f}$ of (9.1); moreover, $\|V^+ \vec{f}\|_2$ is the minimum value of the Euclidean norm $\|\vec{x}\|_2$ over all least-squares solution vectors $\vec{x}$ for such a system.} Combining (9.2) and fact 10.1 implies {\bf Fact 10.2.}  {\it The set of all least-squares solutions to the system of (9.1) coincides with the set of vectors representable as $V^+\vec{f} + \vec{z}$, where $\vec{z} \in {\bf N}$(V). }Now, suppose that we seek a set of least-squares solutions to the linear system $V\vec{x}=\vec{f}$, where $V$ is a singular Toeplitz-like matrix. Then algorithm 8.1 and theorem 8.1 give us ${\bf N}$(V). Due to fact 10.2, it remains to compute $V^+\vec{f}$. In this and next sections, we will use the following definition and auxiliary result. {\bf Definition 10.4.} For a given triple $(A, \vec{v}, m)$ of an $n \times n$ matrix $A$, an $n$-dimensional vector $\vec{v}$, and a natural $m$, the Krylov matrix $K(A,\vec{v}, m)$ is the $n \times m$ matrix with the $i$-th column $A^{i-1}\vec{v}$, $i=1, 2, \ldots, m$. The main algorithm of [P90/92] supports the following result. {\bf Fact 10.3.} {\it For an $n \times n$ Toeplitz-like matrix $A$, for a vector $\vec{v}$, and for $m=O(n)$, the Kryvol matrix $K(A, \vec{v}, m)$ can be computed at the deterministic cost $C_n$ of (1.2).}Seeking $V^+\vec{f}$, we may assume that the matrix $V$ is symmetric since so are the matrices $V^TV$ and $VV^T$ and since we have the next simple result: {\bf Fact 10.4.} $V^+=(V^TV)^+V^T=V^T(VV^T)^+$. An alternative reduction to the symmetric case relies on the following matrix equation: $$ \left( \begin{array}{lr}O & V \\ V^T & O \end{array}  \right)^+ = \left( \begin{array}{lr}O & (V^T)^+ \\ V^+ & O \end{array}  \right) \ . \eqno (10.1) $$ The next fact can be used as a basis for our first algorithm for   the computation of $V^+\vec{f}$, where $V$ is symmetric.{\bf Fact 10.5} [P90]. {\it Let $A$ be an $n \times n$  symmetric matrix, let $C_A(x)=\sum_{i=n-r}^n c_{n-i}x^i$ be its characteristic polynomial, and let $r= {\rm rank} \ A$. Then $$ c_r \neq 0, \  c_rA^+=\sum_{i=1}^{r-1} ((c_{r-1}/c_r)c_i -c_{i+1}) A^{r-i} + (c_{r-1}/c_r)A^r \ .  %\eqno (10.1) $$  }{\bf Algorithm 10.1.} {\bf Input:} a displacement generator of a length $\ell=O(1)$ for an $k \times h$ Toeplitz-like matrix $V$ and an $\ell$-dimensional vector $\vec{f}$. {\bf Output:} the vector $\vec{x}=V^+\vec{f}$. {\bf Computations.} \begin{enumerate} \item  %no 1Compute the characteristic polynomial $C_U(x)= \sum_{i=0}^{2r} c_i(U) x^{k+h-i}$ of the $(k+h) \times (k+h)$ symmetric matrix $$U=\left( \begin{array}{lr} O & V \\V^T & O \end{array} \right) ,  \eqno (10.2) $$where $r = {\rm rank} \  V$. (Note that rank $U=2r$ and $c_{2r}(U)\neq 0$ since $U$ is a symmetric matrix having a rank $2r$.) \item %no 2 Compute the vector $\left( \begin{array}{c} \vec{0}\\ \vec{x}\end{array} \right) = U^+ \left( \begin{array}{c} \vec{f}\\ \vec{0}\end{array} \right)$;  %\item %no 3 Compute and output the vector $\vec{x}$. %=(\frac{1}{c_{2r}(U)})\vec{y}$. \end{enumerate}{\it Correctness} of this algorithm is immediately verified, since $U^+=\left( \begin{array}{lr} O & (V^T)^+ \\V^+ & O \end{array} \right) $ [cf. (10.1)] and, consequently, $U^+ \left( \begin{array}{c} \vec{f}\\ \vec{0}\end{array} \right)=  \left( \begin{array}{c} \vec{0}\\ V^+\vec{f}\end{array} \right)$. Fact 10.5 reduces the computation of the vector $\left( \begin{array}{c} \vec{0}\\ \vec{x}\end{array} \right) = U^+\left( \begin{array}{c} \vec{f}\\ \vec{0}\end{array} \right)$ to the computation of the characteristic polynomial $C_U(x)$ and of the Krylov matrix $K(U,\left(\begin{array}{c} \vec{f}\\ \vec{0}\end{array} \right), r+1)$ (compare definition 10.4). By applying the main algorithm of [P90/92] and fact 10.3, we perform both of these computations at the deterministic cost $C_n$ of (1.2) and thus arrive at the following theorem. {\bf Theorem 10.1.} {\it A least-squares solution $\vec{x}=V^+\vec{f}$ to a Toeplitz-like linear system $V\vec{x}=\vec{f}$ of $k$ equations with $h$ unknowns can be computed at deterministic cost bounded by $C_n$ of (1.2) for $n=k+h$.} In the remainder of this section, we will describe a randomized Las Vegas parallel algorithm for computing $V^+$, which can be immediately extended to the computation of $V^+\vec{f}$. We will start with an auxiliary result. {\bf Fact 10.6} (cf. lemma 2.1 of [B-IG74]). {\it Let the columns of a matrix $B$ form the basis for the null-space ${\bf N}$(A) of a matrix $A$ and let the columns of a matrix $C$form a basis for ${\bf N}$$({\rm A}^T)$. Let $$W= \left( \begin{array}{lr}A & C \\ B^T & O \end{array}  \right) \ . \eqno (10.2) $$ Then $W$ is a nonsingular matrix, and $$ W^{-1} = \left( \begin{array}{lr}A^+ & (B^T)^+ \\ C^+ & O \end{array}  \right) \ .  \eqno (10.3)$$} Fact 10.6 enables us to compute a displacement generator of $V^+$ as follows. {\bf Algorithm 10.2.} {\bf Input:} a field ${\bf F}$, a finite set $S$ in ${\bf F}$ or in its extension field, and a displacement generator of a length $\ell$ for a Toeplitz-like matrix $V$. {\bf Output:} a displacement generator of a length at most $5 \ell + 18$ for $V^+$. {\bf Computations:} \begin{enumerate}\item %no 1 Apply algorithm 8.1 twice, for $A=V$ and $A=V^T$, in order to compute displacement generators of lengths at most $ 2 \ell + 10$ for two matrices $B$ and $C$, whose columns form bases for the null-spaces ${\bf N}$(V) and ${\bf N}$ $({\rm V}^T)$, respectively.\item %no 2 Compute a displacement generator of a length at most $5 \ell +16$ for the matrix $W$ of (10.2), where $A=V$.  \item %no 3 Apply the algorithm supporting theorem 6.1 in order to compute a displacement generator $(G^*, H^*)$ of a length $\ell^* \leq 5 \ell +18$ (cf. fact 2.4) for the matrix $W^{-1}$ of (10.3), where $A=V$. Delete the last rows of the matrices $G^*$ and $H^*$ so as to arrive at $n \times \ell^*$ matrices $G$ and $H$. Output the pair $G$, $H$ as a displacement generator of length $\ell^*$ for $V^+$.\end{enumerate}{\it Correctness} of algorithm 10.2 follows from facts 7.3 and 10.6. %$\left( \begin{array}{lr} %E & 0\\ 0 & 0 \end{array} \right)$, $\left( \begin{array}{lr} %0 & E\\ 0 & 0 \end{array} \right)$, and $\left( \begin{array}{lr} %0 & 0\\ E & 0 \end{array} \right)$  The randomized {\it computational cost} of performing algorithm 10.2 is bounded by $C_n$ of (1.2), since stage 2 is performed at a lower cost whereas stages 1 and 3 are performed at the computational cost estimated in theorems 8.1 and 6.1, respectively, in the case of $p=0$. The bounds of remark 5.4 on the number of random parameters and on the probability of failure apply.  {\bf Theorem 10.2.} {\it At a randomized Las Vegas computational cost bounded by $C_n$ of (1.2), one may compute the Moore-Penrose generalized (or pseudo) inverse of a $k \times h$ real Toeplitz-like matrix (for $k+h=n$) given with its displacement generator of a length $\ell$.} Theorem 10.2 immediately gives us an alternative (although this time randomized) algorithm for computing $V^+\vec{f}$. {\bf Remark 10.1.} We may substantially simplify stage 1 of algorithm 10.2 by shifting from the input matrix $V$ to any of the symmetric matrices $V^TV$, $VV^T$, or $\left( \begin{array}{lr} O & V\\ V^T & O \end{array} \right)$ (cf. (10.1) and fact 10.4). \vspace{.2in} {\bf 11. Extension to computations with other structured matrices and with general matrices.} We have the following remarks on the extension of our results to computations with other classes of matrices. {\bf Remark 11.1.} The presented algorithms and all the asymptotic estimates of theorems 5.1, 6.1, 7.1, 7.2, 8.1, 9.1, 10.1 and 10.2 for parallel arithmetic complexity can be extended to the respective computations with Hankel-like and with Toeplitz-like + Hankel-like matrices. The extension to the Hankel-like matrices immediately follows from their definition (see definition 2.2). The extension to the Toeplitz-like + Hankel-like matrices is obtained by using their representation by the associated short displacement generators (compare remark 2.2 and [KKM79], [BP93], [BP94]). {\bf Remark 11.2.} It was shown in [P89/90] (compare section 12 of chapter 2 of [BP94]) that pre- and/or post-multiplication by appropriate (nonsingular) Vandermonde matrices transforms the associated short (displacement/scaling) generators for structured matrices of the classes of Cauchy-like (Hilbert-like) and Vandermonde-like matrices into short displacement generators of some Toeplitz-like matrices (compare our remark 2.2). Since multiplication of a Vandermonde matrix and of its inverse by a vector reduces to multipoint evaluation and interpolation for the associated univariate polynomial, these operations can be performed at a lower parallel cost over any field of constants (compare [BP94]). This observation enables us to extend the results of theorems 5.1, 6.1, 7.1, 7.2, 8.1 and 9.1 to the cases of Cauchy-like and Vandermonde-like input matrices. {\bf Remark 11.3.} The techniques of [W86], [KP91], [KP92], and [BP94] reduce several fundamental parallel computations with general matrices to ones with Toeplitz, Toeplitz-like, and/or Toeplitz-like + Hankel-like matrices. (This includes solution of nonsingular  and singular linear systems of equations and computation of the inverse, the determinant, the rank, the null-space, the minimum polynomial, a nonsingular submatrix of the maximum size, and various factorizations for a general input matrix.)  Consequently, the present paper should affect the known record estimates for the latter computational problems. The results of the work on this subject  (continuing the study in [KP91], [KP92], and [KS91]) are planned to be published in [KP,a] and/or [BP,a]. \vspace{.3in} {\bf 12. Impact on the computation of Pad\'{e} table, linear recurrence coefficients (Berlekamp-Massey problem), polynomial gcd, lcm, and extended Euclidiean scheme.} Several computational problems of practical importance are reduced to computations with structured matrices. Particularly major are some problems reducible to computing the rank of an $N \times N$ Toeplitz-like matrix and to  $ \hbox{LIN} \cdot \hbox{SOLVE} \cdot \hbox{T} (n)$,  for appropriate $N$ and $n$. In particular, due to facts 2.3 and 2.10, we have the following result. {\bf Theorem 12.1.} {\it The computation of a $(k,h)$ entry of the Pad\'{e} table for a given formal  power series can be performed over any field ${\bf F}$ of constants, having any characteristic $p$, at deterministic parallel cost bounded by the sum of the cost bound $O(\log k, k \log \log k)$ and $C_h$ (that is, of $C_n$ of (1.2) for $n=h$), if $p=0$, and at randomized Las Vegas parallel cost being the sum of the deterministic bound $O(\log k, k \log \log k)$ with the randomized complexity estimates of theorems 5.1, 7.1 and 7.2  for $\ell = 2$ and $n=h$, if $p>0$.} A number of well-known computational problems is equivalent or is reduced easily to computation of a fixed entry of a Pad\'{e} table. In particular, this includes the computation of the coefficients of a linear recurrence of a length at most $n$ from its first $2n$ terms (Berlekamp-Massey problem), compare [BP94], page 48.We will next introduce some additional definitions and auxiliary results (compare [AHU74], [G84], [BP94]), which we will use in order to show some further extensions. \vskip2mm{\bf Definition 12.1} For a pair of nonzero polynomials $u(x)=\sum_{i=0}^n u_i x^i$ and  $v(x)=\sum_{i=0}^m v_i x^i$, their {\it greatest common divisor (gcd)}, $g(x)=gcd(u,v)=\sum_{i=0}^{\delta} g_i x^i$, $g_{\delta}=1$, is the monic common divisor of the polynomials $u(x)$ and $v(x)$ having the largest degree $\delta$. The least common multiple (lcm) of $u(x)$ and $v(x)$ is $u(x)v(x)/g(x)$. Extended Euclidean scheme consists of 4 sequences $\{q_i(x)\}, \  \{r_i(x)\}, \ \{s_i(x)\}$, and $\{t_i(x)\}$, $i=0, 1, \ldots, \ell +1$, where $s_0(x)=t_1(x)=1$, $s_1(x)=t_0(x)=0$, $r_0(x)=u(x)$, $r_1(x)=v(x)$, $$ r_i(x)=r_{i-2}(x)-q_{i-1}(x)r_{i-1}(x) \ , \eqno (12.1) $$$$ s_i(x)=s_{i-2}(x)-q_{i-1}(x)s_{i-1}(x) \ , \eqno (12.2) $$$$ t_i(x)=t_{i-2}(x)-q_{i-1}(x)t_{i-1}(x) \ , \eqno (12.3) $$$$\deg r_i(x) < \deg r_{i-1}(x) \ , \eqno (12.4)$$ for $i=2, 3, \ldots, \ell$, $r_{\ell +1} \equiv 0$, and consequently, $r_{\ell}(x)=\ell .c.(r_{\ell})gcd(u(x), v(x)) =\ell .c.(r_{\ell})g(x)$, $\ell .c.(r_{\ell})$ denoting the leading coefficients of $r_{\ell}(x)$. From (12.1)--(12.4), one may obtain the following result by induction on $i$ (compare [G84], [BP94], page 40).  \vskip2mm {\bf Fact 12.1.} {\it The following relations hold: $$s_i(x)r_0(x) + t_i(x)r_1(x) = r_i(x), \  \  i=1, \ldots, \ell\ , \eqno (12.5) $$$$ \deg r_{i+1}(x) < \deg r_i(x), \ \  i=1, \ldots, \ell -1 \ , $$$$ \deg t_{i+1}(x) > \deg t_i(x) = \deg (q_{i-1}(x) t_{i-1}(x)) = \deg u(x) - \deg r_{i-1}(x), \  \  i=2, 3, \ldots, \ell \ , $$ where $\deg s_{i+1}(x) = \deg t_{i+1}(x)= \infty$.} {\bf Definition 12.2.} The polynomials $s(x)=s_{\ell}(x)$ and $t(x)=t_{\ell}(x)$ are called the {\it cofactors} of $u(x)=r_0(x)$ and $v(x)=r_1(x)$; the triple $\{ g(x), \  s(x), \  t(x)\}$, of the gcd and the cofactors, is called the {\it extended gcd} of $u(x)$ and $v(x)$. {\bf Fact 12.2} [G84]. {\it \begin{description}\item[a)] Let $n>m$. Then $k=\deg r_i(x)$ for some integer $i$, $1\leq i \leq \ell$, if and only if the matrix $R_k$ of definition 2.9 is nonsingular. \item[b)] Let $k=\deg r_i(x)$, let $c_i$ denote the leading coefficient of $r_i(x)$, let $s_i(x) = c_i\sum_{j=0}^{n-k} y_jx^j$, let    $t_i(x) = c_i\sum_{j=0}^{m-k} z_jx^j$, and let $\vec{y}=\vec{y}(i)$ and $\vec{z}=\vec{z}(i)$ denote the coefficient vectors of the polynomials $c_i^{-1}s_i(x)$ and $c_i^{-1}t_i(x)$, respectively. Then $R_k \left[ \begin{array}{c}\vec{y}\\ \vec{z} \end{array} \right]= \left[ \begin{array}{c} 0  \\ \vdots \\ 0 \\ 1 \end{array}\right]$.\end{description} } Due to facts 12.1 and 12.2, the computation of the gcd and its cofactors (the extended gcd problem) is reduced to computing the rank $r$ of a Sylvester matrix and solving a nonsingular subresultant  linear system of $r$ equations  (see our appendix A on an alternative reduction of the gcd computation to the computation of an entry of a Pad\'{e} table). Due to fact 2.11, Sylvester and subresultant matrices have their displacement generators of lengths at most 2 readily available, so that we may again apply our results of sections 5--7. The computation of the gcd is immediately extended to the computation of the least common multiple (lcm), (see definition 12.1). We arrive at the following results. {\bf Theorem 12.2.} {\it The greatest common divisor (gcd) and the least common multiple (lcm) of 2 polynomials of degrees at most $n$ as well as their extended gcd can be computed over any field of constants at the parallel cost speficied in theorem 9.1. The same result applies to the Berlekamp-Massey problem of the recovery of the coefficients of a linear recurrence sequence of length at most $n$ from its first $2n$ terms. The resultant, $\det R$, of 2 polynomials of degrees $k$ and $h$, respectively, for $k+h \leq n$, can be computed at the cost bounded according to theorem 6.1.} Similarly to the computation of the extended gcd, we may proceed with the computation of the coefficient vectors $\vec{y}=\vec{y}(i)$ of $c_i^{-1}s_i(x)$ and $\vec{z}=\vec{z}(i)$ of $c_i^{-1}t_i(x)$ (compare fact 12.3) and, consequently, of the coefficients of $c_i^{-1}r_i(x)$ [see (12.5)], for any $i$. If the degrees of $r_i(x)$ are known, we may compute all entries of the extended Euclidean scheme (12.1)--(12.4) by applying fact 12.2 and the construction of [P89], page 1478, which we will briefly sketch now. We will assume, for simplicity, that $n\geq m$, $m+1$ is an integer power of 2, $m=2^h-1$, and that all the quotient polynomials $q_i(x)$ are linear. Then, the computation of all the entries of the extended Euclidean scheme is reduced to the inversion of the subresultant matrices $R_i$ for $i$, $i=0, 1, \ldots, 2^h-1$, and to simpler operations of multiplication of pairs of polynomials and of Toeplitz-like matrices by vectors. We perform the inversion in $h$ stages. At the $j$-th stage, we invert the matrices $R_k$ for $k=(2s-1)2^{h-j-1}$, $s=1, 2, \ldots, 2^j$; $j=0, 1, \ldots, h$.  At the 0-th stage, we simply apply theorem 6.1; at other stages we first apply the Sherman-Morrison-Woodbury formula (GL89], page 51), which enables us to use the matrix inverses computed at the previous stages. Thus, at the $j$-th stage for $j>0$, we deal with computing short displacement generators for at most $2^j$ products of pairs of Toeplitz-like matrices having sizes less than $(m+n) \times (m+n)$ and with inversion of nonsingular Toeplitz-like matrices of sizes $\ell_j \times \ell_j$, $$\ell_j=(m+n)/2^{j+1} \ .  \eqno (12.6) $$ We apply fact 2.4 and theorem 6.1 in order to perform these operations and to estimate their computational cost. By summarising the resulting estimates and by using Brent's principle, we will extend the estimates of theorem 12.1 to bound the complexity of computing the extended Euclidean scheme. In this extension, the parallel time bounds of theorem 12.1 shall increase by factor $\log n$, whereas the processor bounds shall decrease by the same factor.  An alternative algorithm of [BG90] leads to the same estimates for the complexity of computing the extended Euclidean scheme, provided again that the degrees of all the remainders $r_i(x)$ are known. (The estimates have been stated in [P89] and [BG90] only in the case of the computation over the fields of characteristic $p=0$ and without optimizing the processor bounds, but these estimates can be easily extended and refined.)  The algorithm of [BG92] (reproduced as algorithm 2.10.1 on pages 170--171 of [BP94]) enables us to remove the latter restriction that the degrees are known. The algorithm of [BG92] is reduced to block triangular factorization of a Hankel-like matrix and is quite involved. It can be implemented  over any field of constants; in the case of numerical computations, it enables good control of computational precision. The algorithm consists  of $K=O(\log n)$ stages (compare the proofs of proposition 4.5.3 and corollary 4.5.2 on pages 324 and 325 of [BP94]). The computational cost of performing the $j$-th stage, is dominated by operations with $O(2^j)$ Hankel or Hankel-like matrices of size $\ell_j \times \ell_j$, $j=0, \ldots, K-1$ [cf. (12.6)]. These operations include  computing the rank of a Hankel matrix and/or the inverse of a nonsingular Hankel matrix, as well as a pair of polynomials generating a fixed Hankel matrix (see proposition 2.9.1 on page 158 of [BP94]). The latter operation is reduced to computing the $(k_j,h_j)$ entry of a Pad\'{e} table for $k_j +h_j \leq \ell_j$. By combining theorems 6.1, 7.1, 7.2, and 12.1, we arrive at the desired extension of the estimates of theorem 12.1. The next theorem summarizes the resulting complexity estimates (whose further elaboration is left for the reader). {\bf Theorem 12.3.} {\it For 2 polynomials $u(x)$ and $v(x)$ of degrees at most $n$, all the entries $r_i(x)$, $q_i(x)$, $s_i(x)$ and $t_i(x)$ of the associated extended Euclidean scheme (12.1)--(12.4) can be computed over any field ${\bf F}$ of constants, having a   characteristic $p$, at the deterministic cost $O((\log n)^3, (n^2 \log \log n)/(\log n)^2)$ if $p=0$ and [with a probability of failure at most $2^{-\beta}$ for any positive $\beta=O(\log n)$] at a Las Vegas randomized cost, bounded by $O((\log n)^3, (n^2 \log \log n)/((\log n)\log \min \{ |{\bf F}|, n\}))$ if $p>n$ and by \\$O((\log n)^3 (\gamma(n,p))^2, (n^2 \log \log n)/((\gamma(n,p))^2(\log n)\log \min \{ |{\bf F}|, n\}))$ if $0< p <n$, where $\gamma(n,p)=\lceil \log n/\log p \rceil$. The randomized computations involve $n^{O(1)}$ random parameterschosen from a fixed finite set $S$ (in ${\bf F}$ or in its algebraic extension field) that has a cardinality $|S|=n^{O(1)}$.} \vskip2mm\noindent{\bf Acknowledgements.}  This work has been supported by NSF grant CCR 9020690 and PSC CUNY Awards 665301 and 666327. Erich Kaltofen pointed out a flaw in my firstproof of theorem 5.1. Igor Shparlinski furnished a reference on application of our results to algebraic codes. Two referees provided many helpful comments and suggested to expand the bibliography and the material on the historical background, as well as on applications and extensions of our algorithm 5.1.\vskip4mm\begin{center}{\bf REFERENCES}\end{center}\begin{description}\item [AHU74] A.V. Aho, J.E. Hopcroft, J.D. Ullman, {\it The Designand Analysis of Computer Algorithms}, Addison-Wesley, Reading,Massachusetts, 1974. \item[AGr88]  G.S. Ammar, W.G. Gragg,  Superfast Solution of Real Positive  Definite Toeplitz Systems,  {\it  SIAM J.  Matrix Anal.  Appl.,} 9, 1,  61--76,  1988. \item [BA80]  R.R. Bitmead, B.D.O. Anderson,  Asymptotically Fast Solution of Toeplitz and Related Systems of Linear Equations,  {\it  Linear Algebra Appl.}, 34, 103--116,  1980.\item[Be68]		E. R. Berlekamp, {\it Algebraic Coding Theory,} McGraw-Hill, New York, 1968.\item [Be84]  S. Berkowitz,  On Computing the Determinant in Small Parallel Time Using a Small Number of Processors,  {\it  Inform. Proc. Lett.,} 18, 3,  147--150,  1984.\item[BG90]  D. Bini, L. Gemignani, On the Euclidean Scheme for Polynomials Having Interlaced Real Zeros,    {\it Proc.  2nd Ann. ACM Symp. on ParallelAlgorithms and Architectures},  254--258, ACM Press, New York, 1990. %\pndt % p904\item[BG92]  D. Bini, L. Gemignani,  Fast Parallel Computation of the Remainder Sequence via Bezout and Hankel Matrices,  {\it Quaderno 79, Dipartimento di Matematica, Universit\`a di Parma,} Parma, Italy, 1992 (to appear in {\it SIAM J. Comput.}). %\pndt % p879\item [BGH82]  A. Borodin, J. von zur Gathen, J. Hopcroft,  Fast Parallel Matrix and GCD Computation,  {\it  Information and Control}, 52, 3, 241--256,  1982.  \item [BGY80] R. P. Brent, F. G. Gustavson, D. Y. Y.Yun, Fast Solution of Toeplitz Systems of Equations  and Computation  of  Pad\'{e} Approximations, {\it J. of Algorithms}, 1, 259-295, 1980.\item[B-IG74] A. Ben-Israel, T. N. E. Grevill, {\it Generalized Inverses: Theory and Applications}, John Wiley, New York, NY, 1974.  \item[BK87]		A. Bruckstein, T. Kailath, An Inverse Scattering Framework for Several Problems in Signal Processing, {\it IEEE ASSP Magazine,} 6-20, January 1987.\item[B-OT90] M. Ben-Or, P. Tiwari, Simple Algorithm for  Approximating All Roots of a Polynomial with Real Roots, {\it J. of Complexity}, 6, 4, 417--442, 1990.\item[BP91]  D. Bini, V. Y. Pan,  Parallel Complexity of Tridiagonal Symmetric Eigenvalue Problem, {\it Proc. 2nd Ann. ACM-SIAM Symposium on Discrete Algorithms,} 384--393, ACM Press, New York, and SIAM Publications, Philadelphia, Pennsylvania, 1991.\item [BP92/93]  D. Bini, V. Y. Pan,  Improved Parallel Polynomial Division,  {\it Proc. 33rd Ann. IEEE Symp. on Foundation of Computer Science,} 131--136, IEEE Computer Society Press, 1992, and {\it  SIAM J. Comput.}, 22,  3, 617--626, 1993. \item [BP93]D. Bini, V.Y. Pan,  Improved Parallel Computations with Toeplitz-like and Hankel-like Mtrices, {\it Linear Algebra Appl.}, 188, 189, 3--29, 1993. [Proceedings version, {\it Proc. ACM-SIGAM Int. Symp. on Symbolic and Alg. Comp. (ISSAC '93)}, (M. Bronstein Editor) 193--200, ACM Press, New York, 1993.] \item [BP94]D. Bini, V.Y. Pan, {\it Matrix and Polynomial Computations, Volume 1: Fundamental Algorithms}, Birkhauser,Boston, 1994. \item [BP,a]D. Bini, V.Y. Pan, {\it Matrix and Polynomial Computations, Volume 1: Fundamental Algorithms} ($2^{nd}$ edition), to appear.\item [Bre74]R.P. Brent, The Parallel Evaluation of General Arithmetic Expressions, {\it J. of ACM}, {\bf 21},2, 201-206, 1974.  \item[Bro71] W. S. Brown, On Euclid's Algorithm and the Computation of Polynomial Greatest Common Divisors, {\it J. of ACM}, 18, 478--504, 1971.\item[BT71]  W.S. Brown, J.F. Traub,  On Euclid's Algorithm and the Theory of Subresultants,  {\it  J. of ACM,} 18 ,  4 ,  505--514, 1971. \item  [Ch85]  A.L. Chistov,  Fast Parallel Calculation of the Rank of Matrices over a Field of Arbitrary Characteristics,  { \it  Proc. FCT 85, Lecture Notes in Computer Science}, 199, 63--69, Springer, Berlin, 1985.  \item[ChK91]  J. Chun, T. Kailath,  Divide-and-Conquer Solution of Least-Squares Problems for Matrices with Displacement Structure,   {\it SIAM J.   Matrix Anal.  Appl.,} 12, 1, 128--145, 1991.\item [CK91]  D. G. Cantor,  E. Kaltofen, On Fast Multiplication ofPolynomials over Arbitrary Rings, {\it Acta Inf.}, 28, 7, 697--701, 1991. \item[CKL-A87]  J. Chun, T. Kailath, H. Lev-Ari,  Fast Parallel Algorithm for QR-factorization of Structured Matrices, {\it  SIAM J.   Sci.  Stat.Comput.,}  8, 6,  899--913,  1987.\item[Col67] G. E. Collins, Subresultants and Reduced Polynomial Sequences, {\it J. of ACM}, 14, 128--142, 1967. \item[Cs76]  L. Csanky,  Fast Parallel Matrix Inversion Algorithm,  {\it  SIAM J. Comput.,}   5, 4,  618--623,  1976.\item[CW90]D. Coppersmith, S. Winograd, Matrix Multiplication via ArithmeticProgressions, {\it J. of Symbolic Computation}, 9, 3, pp. 251-280,1990.\item[dH87]  F.R. deHoog,  On the Solution of Toeplitz Systems,  {\it  Linear Algebra  Appl.,}  88/89,  123--138,  1987.  \item [DL78] R.A. Demillo, R.J. Lipton, A Probabilistic Remark on AlgebraicProgram Testing, {\it Information Process. Letters}, 7, 4, 193--195, 1978. \item[EG88]  D. Eppstein, Z. Galil,  Parallel Algorithmic Techniques for Combinatorial Computation,  {\it  Annual Review of Computer Science,} 3,  233--283,  1988.\item [G84]  J. von zur Gathen,  Parallel Algorithms for Algebraic Problems,{\it  SIAM J.  Comput.},  13, 4,  802--824,  1984.   \item[GG94] M. H. Gutknecht, W. B. Gragg, Stable Look-ahead Versions of the Euclidean and Chebyshev Algorithms, IPS Research Report 94-04, ETH-Zentrum, Z\"urich, Switzerland, 1994. \item[GP89]  Z. Galil, V. Y. Pan,  Parallel Evaluation of the Determinant and of the Inverse of a Matrix,  {\it  Inf. Proc. Lett., } 30,  41--45,  1989. \item[Gr72] W. B. Gragg, The Pad\'{e} Table and Its Relation to Certain Algorithms of Numerical Analysis, {\it SIAM Review}, 14, 1, 1--62, 1972. \item[J92] J. Ja Ja, {\it An Introduction to Parallel Algorithms,} Addison-Wesley, Mas\-sa\-chu\-setts, 1992.  \item[K74]			T. Kailath, A View of Three Decades of Linear Filtering Theory, {IEEE Trans. Inform. Theory,} 20, 146--181, 1974.\item[K87]  T. Kailath,  Signal Processing Applications of Some Moment Problems, {\it  Proc. AMS Symp. in Applied Math.,} 37,  71--100,  1987.  \item [K94] E. Kaltofen, Asymptotically Fast Solution of Toeplitz-like Singular Linear Systems, {\it Proc. ACM-SIGSAM Internat. Symp. on Symbolic and Algebraic Comput. (ISSAC' 94)}, edited by J. von zur Gathen and M. Giesbrecht;  297--304, ACM Press, New York, NY, 1994.\item [K,a] E. Kaltofen, Analysis of Coppersmith's Block Wiedemann Algorithm for the Parallel Solution of Sparse Linear Systems, {\it Math. Comput.}, to appear.\item [Ka25] S. Kakeya, On Fundamental Systems of Symmetric Functions, {\it Japan J. of Math.}, 2, 69--80, 1925. \item[Ka27] S. Kakeya, On Fundamental Systems of Symmetric Functions, {\it Japan J. of Math.}, 4, 77--85, 1927.\item[KC89]  T. Kailath, J. Chun, Generalized Gohberg-smencul Formulas for Matrix Inversion, {\it Operator Theory: Advances and Apllications}, 40, 231--246, 1989.\item[KKM79]  T. Kailath, S.-Y. Kung, M. Morf,  Displacement Ranks of Matrices and Linear Equations,  {\it  J. Math. Anal. Appl.},  68, 2,  395--407,  1979. \item [KP91]E. Kaltofen, V.Y. Pan, Processor Efficient Parallel Solution of LinearSystems over an Abstract Field, {\it Proc. 3rd Ann. ACM Symp. on ParallelAlgorithms and Architectures}, 180-191, ACM Press, New York, 1991.\item [KP92]E. Kaltofen, V.Y. Pan, Processor Efficient Parallel Solutionof Linear Systems II. The Positive Characteristic and Singular Cases,{\it Proc. 33rd Ann. IEEE Symp. on Foundation of Computer Science}, 714-723,IEEE Computer Society Press, 1992.\item [KP94]E. Kaltofen, V.Y. Pan,Parallel Solution of Toeplitz and Toeplitz-like Linear Systems over Fieldsof Small Positive Characteristic,{\it Proc. First Intern. Symposium on ParallelSymbolic Computation} (PASCO '94), Linz (edited by Hoon Hong),Austria (Sept. 1994),{\it Lecture Notes in Computing,} {\bf 5}, 225-233,World Scientific Publishing Company, Singapore,1994.\item[KP,a] E. Kaltofen, V. Y. Pan, to appear. \item [KR90]R. Karp, V. Ramachandran, A Survey of Parallel Algorithms for SharedMemory Machines, {\it Handbook for Theoretical Computer Science} (J. vanLeeuwen Editor), 869-941, North Holland, Amsterdam, 1990.\item[KS91] E. Kaltofen, B. D. Saunders, On Wiedemann's Method for Solving Sparse Linear Systems,  {\it Proc. AAECC-5,   Lecture Note in Computer Science}, 536, 29--38, Springer, Berlin, 1991.  \item[KVM78]		T. Kailath, A. Viera, M. Morf, Inverses of Toeplitz Operators, Innovations, and Orthogonal Polynomials, {\it SIAM Rev.,} 20, 1, 106--119, 1978.\item[L-AK84]		H. Lev-Ari, T. Kailath, Lattice Filter Parametrization and Modelling of Nonstationary Processes, {\it IEEE Trans. on Information Theory,} IT-30, 1,  2--16, 1984.\item[L-AKC84]		H. Lev-Ari, T. Kailath, J. Cioffi, Least Squares Adaptive Lattice and Transversal Filters; a Unified Geometrical Theory, {\it IEEE Trans. on Information Theory,} IT-30, 222--236, 1984.\item[Ma75]		J. Makhoul, Linear Prediciton: A Tutorial Review, {\it Proc. IEEE,} 63, 4, 561--580, 1975.\item [Morf80]  M. Morf,  Doubling Algorithms for Toeplitz and Related Equations, {\it  Proc. IEEE Internat. Conf. on ASSP,} 954--959, IEEE Computer Society Press, 1980. \item[Mus81] B.R. Musicus, Levinson and Fast Choleski Algorithms for Toeplitzand Almost Toeplitz Matrices, \rm Internal Report, {\it Lab. of Electronics, M.I.T.,} Cambridge, Massachusetts,  1981.\item[N27] K. Nakamura, On the Representation of Symmetric Functions by Power-Sums Which Form the Fundamental System, {\it Japan J. Math.}, 4, 87--92, 1927.  \item [P85]V.Y. Pan, Fast and Efficient Parallel Algorithms for the Exact Inversionof Integer Matrices, {\it Proc. 5th Ann. Conference on Foundation ofSoftware Technology and Theoretical Computer Science, Lecture Notes inComputer Science}, 206, 504-521, Springer, Berlin, 1985.\item [P87]V.Y. Pan, Complexity of Parallel Matrix Computations, {\it Theoretical Computer Science}, 54, 65--85, 1987. \item[P88]  V. Y. Pan,  New Effective Methods for Computations with Structured Matrices, {\it  Technical Report 88-28,  Computer Science Dept.,  SUNY Albany,} Albany, New York,  1988. \item[P89]  V. Y. Pan,  Fast and Efficient Parallel Evaluation of the Zeros of a Polynomial Having Only Real Zeros,  {\it  Computers \& Mathematics (with Applications),}  17, 11,   1475--1480,  1989. \item[P89/90] V. Y. Pan, On Computations with Dense Structured Matrices, {\it Proc. ACM-SIGSAM Intern. Symp. on Symbolic and Algebraic Comp.} (ISSAC' 89), 34--42, ACM Press, New York,  1989, and {\it  Math. Comp.},  55, 191,  179--190, 1990.\item[P91]			V. Y. Pan, Complexity of Algorithms for Linear Systems of Equations, in {\it Computer Algorithms for Solving Linear Algebraic Equations (The State of the Art)}, edited by E. Spedicato, NATO ASI Series, Series F: Computer and Systems Sciences, 77, 27-56, Springer, Berlin (1991).\item [P90/92]V. Y. Pan, Parametrization of Newton's Iteration for Computations withStructured Matrices and Applications, {\it Computers \& Mathematics (withApplications)}, {\bf 24}, 3, 61-75, 1992.\item [P92a] V. Y. Pan, Complexity  of Computations with Matrices and Polynomials, {\it SIAM Rev.}, 34, 9, 225--262, 1992.\item [P92c] V. Y. Pan,Parallel Solution of Toeplitzlike Linear Systems, {\it J. of Complexity}, 8, 1--21, 1992.\item[P93]  V. Y. Pan, Decreasing the Displacement Rank of a Matrix,{\it SIAM J.  Matrix Anal. Appl.,} 14, 1, 118--121, 1993.\item[P93a] V. Y. Pan, On Randomized Parallel Computations with General and Toeplitz-like Matrices, manu\-script, 1993. \item [P93b] V. Y. Pan, Concurrent Iterative Algorithm for Toeplitz-like Linear Systems,  {\it IEEE Trans. on Parallel and Distributed Systems}, 4, 5, 592--600, 1993. \item[P95]V.Y. Pan, Optimal and Nearly Optimal Algorithms for Approximating Complex Polynomial Zeros, Computer and Math. (with Applics),  {\em Proc. 27th Ann. ACM Symp. on Theory of Computing}, ACM Press, New York, 1995.\item[P,a] V. Y. Pan, Effective Parallel Computations with Toeplitz and Toeplitz-like Matrices Filled with Integers, to appear in {\it Proc. AMS-SIAM Workshop on Mathematics of Numerical Analysis: Real Number Algorithms} (M. Shub, S. Smale, J. Renegar editors, Park City, Utah, June--July 1995), Amer. Math. Soc., Providence, RI, 1996.\item[P,b] V. Y. Pan, Pad\'{e} Approximation Approach to the Solution of the Key-Equation of the Algebraic Coding Theory, preprint, 1996.  \item [PP95]V. Y. Pan, F. P. Preparata, Work-Preserving Speed-up of Parallel Matrix Computations, {\em SIAM J. on Computing}, 24, 4, 811--821, 1995. \item[PR91] V. Y. Pan, J. Reif, The Parallel Computation of the Minimum Cost Paths in Graphs by Stream Contraction, {\it Information Processing Letters}, 40, 79--83, 1991.   \item[PSa78]  F. P. Preparata, D. V. Sarwate, An Improved Parallel Processor Bound in Fast Matrix Inversion, {\it Inform. Proc. Letters}, 7, 3, 148--149, 1978.\item[PSLT93]V.Y. Pan, A. Sadikou, E. Landowne, O. Tiga, A New Approach to Fast Polynomial Interpolation and Multipoint Evaluation, {\em Computers and Math. (with Applications)}, 25, 9, pp. 25--30, 1993.\item[Q94] M.J. Quinn, {\it Parallel Computing: Theory and Practice,} McGraw-Hill, New York, 1994. \item[R93]  J. H. Reif, An $O(n\log^3 n)$ Algorithm for the Real Root Problem, {\it Proc. of 34rd Ann. IEEE Symp. on Foundations of Computer Science,} 626--635, IEEE Computer Society Press, 1993.\item[R95]			J. H. Reif, Work Efficient Parallel Solution of Toeplitz Systems and Polynomial GCD, {\it Proc. 27th Ann. ACM Symp. on Theory of Computing,} 751--761, ACM Press, New York, NY, 1995.\item[S80]  J.T. Schwartz,  Fast Probabilistic Algorithms for Verification of Polynomial Identities,  {\it  J. of ACM}, 27, 4,   701--717,  1980. \item [Sc82]A. Sch\"{o}nhage, The Fundamental Theorem of Algebra in Terms ofComputational Complexity, manuscript, {\it Dept. of Math., University ofT\"{u}bingen}, T\"{u}bingen, Germany, 1982.\item [Sc93]A. Sch\"{o}nhage, Fast Parallel Computation of Characteristic Polynomialsby Leverrier's Power Sum Method Adapted to Fields of FiniteCharacteristic, {\it Proc. 20th ICALP, Lecture Notes in Computer Science},{\bf 700}, 410-417, Springer, Berlin, 1993.\item[Sh94] V. Shoup, Fast Construction of Irreducible Polynomials over Finite Fields, {\it J. Symb. Computation}, 17, 5, 371--392, 1994. \item [Ste94] 		W. F. Stewart, {\it Introduction to the Numerical Solution of Markov Chains,} Princeton University Press, New Jersey, 1994.\item[V900] K. T. Vahlen, \"{U}ber Fundamentalsysteme f\"{u}r symmetrisch Functionen, {\it Acta Mathematica}, 23, 91--120, 1900. \item [W86] D. H. Wiedemann, Solving Spare Linear Equations Over  Finite Fields, {\it IEEE Trans. on Inf. Theory},  IT-32, 1, 54--62, 1986.  \item[Z79] R.E. Zippel,  Probabilistic Algorithms for Sparse Polynomials,     { \it  Proc.\break EUROSAM 79,   Lecture Notes in Computer Science}, 72,  216--226, Springer, Berlin, 1979.\end{description}\vspace{ .5in}\noindent {\bf Appendix A. Reduction of the polynomial gcd problem to computing a Pad\'{e} approximation.}Let us reduce the computation of the gcd of 2 polynomials   to computing an entry of an associated Pad\'{e} table. With no loss of generality, we may assume that $V(0) \neq 0$.  Then, we apply the following algorithm (cf. e.g. [BP94]): \begin{enumerate} \item %no 1Compute the coefficients $q_0, q_1, \ldots, q_{k+h}$ of the formal power series  $$ Q(x)=U(x)/V(x)=\sum_{i=0}^{\infty} q_ix^i \ . $$\item %no2 Compute the $(k,h)$ entry $(U_-(x), V_-(x))$ of the Pad\'{e} table for $Q(x)$ (see fact 2.10 and theorem 12.1). \item %no3 Compute and output ${\rm gcd} \ (U(x), V(x)) =U(x)/U_-(x)$. \end{enumerate}Correctness of the algorithm is immediately verified. Stages 1 and 3 can be reduced to 2 polynomial divisions.  Stage 2 is reduced essentially to computing the rank of an $h \times h$ Toeplitz matrix, to solving a nonsingular Toeplitz linear system of at most $h$ equations, and to multiplication of a $k \times k$ triangular Toeplitz matrix by a vector, due to fact 2.10. Numerical stability of the computation by this algorithm in the real and complex fields can be improved by using Bezout matrices (see [BG92] or [BP94], section 9 of chapter 2).\vspace{ .5in}\noindent {\bf Appendix B. Decreasing the length of a displacement generator.} Suppose that an $n \times n$ Toeplitz-like matrix $V$ is given with its displacement generator of a length $\ell$ and has  a known displacement rank $\rho < \ell$. Then, at the cost $O(\ell, n\ell)$, one may compute (over any field of constants) a displacement generator of length $\rho$ for $V$ (see [P90/92], proposition A.6). To compute $\rho = {\rm rank} \ (V-ZVZ^T)$, we may apply fact 7.2 to the unit lower triangular Toeplitz matrices $U^T$ and $L$ and then obtain $\rho$ by means of binary search since the largest nonsingular leading principal submatrix of $U(V-ZVZ^T)L$ has size $\rho \times \rho$. Since $\rho < \ell$, the binary search only requires $\lceil \log \ell\rceil$ nonsingularity tests for $i \times i$ matrices, $i < \ell$. If $\ell =O(1)$, the entire computation of the displacement rank $\rho$ of $V$ and of a displacement generator of length $\rho$ for $V$ is performed at a randomized cost $O(\log n, n \log \log n)$. We could have applied such a reduction of the length  of a displacement generator to every Toeplitz-like matrix that we encounter in our computations, but we were able to devise our algorithms without using this procedure in the present paper. The latter procedure, however, is crucial for devising the RNC and work efficient algorithms of [BP94], [R95], and [P,a], proposed for Toeplitz, Toeplitz-like, Pad\'{e}, gcd and lcm computations, which we cited in the introduction. Let us comment on this issue in the remainder of the appendix.All the algorithms [BP94], [R95], [P,a] rely on using order $\log n$ steps of Newton's iteration, repeated both in its algebraic and its numerical versions. For {\it algebraic version}, the displacement rank of the auxiliary matrices is bounded, as proven in [P90/92]. Without extra care, however, the computed displacement generators of some of these matrices generally have lengths as large as $n$, and then the estimated work-complexity of computations with such matrices and, consequently, the overall work-complexity bound grow roughly by factor $n^{\omega -1}$, where $2\leq \omega < 2.376$. To avoid such a growth, one needs to apply either proposition A.6 of [P90/92] or another technique reaching the same goal. This is  done in [BP94] and [P,a], but the entire problem is overlooked in [R95], where the proposed solution relies on a rather naive claim that the $\rho \times \rho$ leading principal submatrix of the matrix $V-ZVZ^T$ is nonsingular for any matrix $V$ having a displacementrank $\rho$ (see [R95], page 753). (This claim is wrong, for instance, for the matrix $\left[ \begin{array}{lr}0 & 0\\  1 & 0 \end{array} \right]$.) Much more serious difficulties arise with performing (order $\log n$) {\it numerical steps} of Newton's iteration in [R95], because such steps involve auxiliary matrices $V$ whose displacement rank generally grows to the level $n$. In this case, proposition A.6 of [P90/92] cannot help, and more advanced and more intricate techniques (using approximations of singular value decompositions of the matrices $V-ZVZ^T$, as well as estimates for the error norms of these approximations) have been applied in [BP94] and [P,a] in order to overcome these difficulties. In contrast to this, J. H. Reif (in [R95]) does not solve such a major problem but simply overlooks it, thus ignoring the resulting growth of the overall work-complexity bound by factor  $n^{\omega -1}$, where $2\leq \omega < 2.376$.\end{document}  