From:	IN%"sad@cunyvms1.gc.cuny.edu" 25-MAR-1996 14:28:51.19To:	IN%"vpan@lcvax.lehman.cuny.edu"CC:	Subj:	This is PARALLEL COMPLEXITY OF TOEPLITZ....INTEGER INPUT AND EXTENSION.Return-path: <sad@cunyvms1.gc.cuny.edu>Received: from vice1.lehman.cuny.edu by lcvax.lehman.cuny.edu (PMDF #12809) id <01I2RBUXRGB4000A75@lcvax.lehman.cuny.edu>; Mon, 25 Mar 1996 14:26 ESTReceived: from CUNYVMS1.GC.CUNY.EDU by vice1.lehman.cuny.edu (5.51/SMI-3.0DEV3) id AA12178; Mon, 25 Mar 96 14:10:42 ESTReceived: by cunyvms1.gc.cuny.edu (MX V3.3 VAX) id 29483; Mon, 25 Mar 1996 14:05:02 ESTDate: Mon, 25 Mar 1996 14:03:16 ESTFrom: sad@cunyvms1.gc.cuny.eduSubject: This is PARALLEL COMPLEXITY OF TOEPLITZ....INTEGER INPUT AND EXTENSION.To: vpan@lcvax.lehman.cuny.eduMessage-id: <0099FDDC.4761B4EC.29483@cunyvms1.gc.cuny.edu>\documentstyle[12pt]{article}%    Set type page length to 50.5pc but adjust %    for the running heads.%%%%%\setlength{\textheight}{50.5pc}%%%%%\addtolength{\textheight}{-\headsep}%    Note: use 5pt instead of \headheight %    here because we want the%    an approximation of the actual visual %    height of the running head,%    assuming eight or ten point font size.%%%%%\addtolength{\textheight}{-5pt}%    Round type page length to an %    even multiple of line spacing (and%    adjust for the first line, which %    gets \topskip instead of%    \baselineskip).%%%%%\normalsize % Ensure proper font size/baselineskip%%%%%\addtolength{\textheight}{-\topskip}%%%%%\addtolength{\textheight}{.5\baselineskip}%%%%%\divide\textheight by \baselineskip%%%%%\multiply\textheight by \baselineskip%%%%%\addtolength{\textheight}{\topskip}%    Set type page width%%%%%\setlength{\textwidth}{30pc}\topmargin-.25in\oddsidemargin.05in \evensidemargin.05in\textheight8.55in  \textwidth6in\newtheorem {algorithm}{\sc \bf Algorithm}[section]\newtheorem {problem}{\sc \bf Problem}[section]\newtheorem {theorem}{\sc \bf Theorem}[section]\newtheorem {proposition}{\sc \bf Proposition}[section]\newtheorem {lemma}{\sc \bf Lemma}[section]\newtheorem {corollary}{\sc \bf Corollary}[section]\newtheorem {fact}{\sc \bf Fact}[section]\newtheorem {remark}{\sc \bf Remark}[section]\newtheorem {definition}{{\sc \bf Definition}}[section]\newtheorem {claim}{\sc \bf Claim}[section]\def\theequation{\arabic{section}.\arabic{equation}}\begin{document}\title{\bf \Large Parallel Complexity ofToeplitz-like Computations with Integer Inputs andExtensions. }\vspace{.2in} \author{ Victor~Y.~Pan \thanks{Supported by NSF Grant CCR 9020690 and PSC CUNY Award No. 666327.} \\ Department of Mathematics and Computer Science\\Lehman College, City University of New York\\Bronx, NY 10468 \\ Internet:  VPAN@LCVAX.LEHMAN.CUNY.EDU}%\thispagestyle{empty}\date{} \pagestyle{myheadings}\markboth{V. Y. Pan}{Parallel Complexity of Polynomial GCD and Toeplitz-like...  }              \maketitle\vspace{.4in} \begin{center} \small\baselineskip=9pt{\bf Abstract\footnote{The results of this paper have been presented at the ACM-SIAM Workshop ''Mathematics of Numerical Analysis: Real Numbers Algorithms'', Park City, Utah, July--August 1995 (cf. [P,c]).}}\end{center} We present parallel algorithms for computing the  rank, the determinant, %$\det A$, and a basis for the null-space of an $n \times n$ Toeplitz or Toeplitz-like matrix $A$ filled with integers, as well as a solution $\vec{x}$ to a linear system $A\vec{x}=\vec{f}$ if the system is consistent. Our algorithms use $O((\log n) \log (n \log \|A\|))$ parallel time and $O(n \log n)$ processors, each capable of performing (in unit time) an arithmetic operation or rounding-off a rational number to the closest integer.  The computations by these algorithms can be performed with the precision of $O(n \log \|A\|)$ bits, which matches the precision required in order to represent the output. The algorithms involve $2n-1$ random parameters for computing $\det~A$ and a single random parameter for other computations. The results are immediately  extended to the evaluation of the gcd, the lcm, and the resultant of two polynomials, to computing Pad\'{e} approximation  and Berlekamp-Massey recurrence coefficients, and to many other related computational problems reducible  to Toeplitz/Toeplitz-like computations and/or to computations with various other structured matrices. The cited processor bounds are less by roughly factor $n$ than ones supported by the known algorithms that run in polylogarithmic arithmetic time and do not use rounding-off to the closest integers. Our algorithm exploits some new techniques as well as our earlier techniques  of variable  diagonal (relating some known algebraic and numerical methods), stream contraction, and the SVD trancation of displacement generators in Toeplitz-like computations; our development and application of these techniques may be of independent interest.\vspace{.2in}\noindent{\bf Key words:} parallel algorithms, Toeplitz matrix computations, Toeplitz-like matrices, polynomial gcd, displacement rank, computational complexity, block Gauss-Jordan decomposition, $p$-adic lifting, Newton-Hensel's lifting, symbolic linear algebra.%\end{abstract}\vspace{.2in}\noindent{\bf 1991 Mathematics Subject Classification:} 68Q22, 68Q25, 68Q40, 65Y20, 47B35, 65F30. \vspace{.23in}%\noindent \begin{center} {\bf 1. Introduction}\end{center}{\bf 1.1.  Toeplitz and Toeplitz-like matrices and some applications.}  %%%  Our goal and background.} The fast version of Euclidean algorithm [AHU74], [BGY80] computes the greatest common divisor (gcd) of 2 polynomials of degrees at most $n$ by using $O(n (\log n)^2)$ arithmetic operations, but to yield substantial parallel acceleration, one has to reduce the problem to the solution of the associated (possibly singular)  Toeplitz-like  (resultant or subresultant) or Toepltz linear system of $N=O(n)$ equations $V\vec{x}=\vec{f}$ [BT71], [BGH82], [G84], [BP94]. The gcd computation is only one (though celebrated) example of numerous major problems of algebraic and numerical computing whose solution is reduced to solving Toeplitz or Toeplitz-like linear systems of equations. The list of such problems includes the evaluation of the  resultant, of the Sturm and subresultant sequences, of the least common multiple (lcm) for a pair of univariate polynomials, computation of the Pad\'{e} approximation of an analytic function, Berlekamp-Massey problem of the recovery of the coefficients of a linear recurrence sequence from its terms, 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], [BT71], [K74], [Ma75], [KVM78], [BGY80], [L-AK84], [L-AKC84], [BK87], [K87], [P91], [P92a], [KP91], [PSLT93], [KP92], [BP94], [Ste94], [P95], [P96,a]. Furthermore, the general reduction techniques of [P90], [He94], [GKO94] enable us to extend the algorithms available  for Toeplitz and Toeplitz-like computations to computations with some other major classes of structured matrices, such as Cauchy-like and Vandermonde-like matrices, also highly important in many areas of computing.Parallel solution of Toeplitz and Toeplitz-like linear systems will be our main topic. For reader's convenience, we will recall the definitions and some basic properties of Toeplitz and Toeplitz-like matrices, at first , now and then, with more details, in section 13. $V=(v_{i,j})$ is an $N \times N$ {\it Toeplitz matrix} if $$ v_{i,j}=v_{i+1,j+1},~~{\rm for}~~i,j=0,1,\ldots, N-2, \eqno (1.1) $$ that is, if the entries of $V$ are invariant in their shifts into the diagonal direction. Toeplitz matrices are easy to store, since such an $N \times N$ matrix is fully represented by the $2N-1$ entries of its first column (or row) and of its last  column (or row). Multiplication of an $N \times N$ Toeplitz matrix by a vector can be reduced to polynomial multiplication modulo $x^{2N-1}$ (see [BP94], p. 133), and then, or more directly, to FFTs, and can be performed by using $O(N\log N)$ arithmetic operations, hereafter referred to as {\it ops}. (Note that for a general matrix, we would need $N^2$ storage and $2N^2-N$ ops, respectively.) Due to the structural properties of Toeplitz matrices, one may solve a nonsingular Toeplitz linear system of $N$ equations by using $O(N(\log N)^2)$ ops [BGY80], [BA80], [Mu81], [dH87], [AGr88], [K94]. The above properties are extended to the class of $N \times N$ {\it Toeplitz-like matrices}, that is, ones representable in the form $$ T=\sum_{i=1}^{\ell} L_iU_i,  \eqno (1.2) $$where $L_i$ and $U_i^T$ are $N \times N$ lower triangular Toeplitz matrices, $U_i^T$ is the transpose of $U_i$, and $\ell$ is a fixed constant, $\ell =O(1)$. Then, we may store $2\ell N$ entries of the first columns of $L_i$ and $U_i$ for $i=1, \ldots, \ell$ in order to represent $T$. Furthermore, we may multiply $T$ by a vector by using $O(\ell N \log N)$ ops, and solve a linear system $T\vec{x} =\vec{f}$ in $O(\ell^2 N (\log N)^2)$ ops if $T$ is nonsingular.We have $\ell < 2$ in (1.2) for Toeplitz matrices, their inverses, resultant and subresultant matrices, $\ell \leq g+h$ for $g \times h$  block matrices with Toeplitz blocks. Furthermore, transposition of a matrix leaves $\ell$ invariant, whereas $\ell$ grows only slowly in multiplicaton and addition/subtration of 2 Toeplitz-like matrices, as well as in their inversion if they are nonsingular.{\bf 1.2. NC and  RNC solutions (some background).} Due to their reduction to Toeplitz/Toeplitz-like linear systems, the several computational problems listed in the previous section, including the gcd, lcm, and resultant computation and Pad\'{e} approximation, can be solved by using $O(N(\log N)^{d_1})$ ops, where $N$ is the input size, and $d_1 \leq 3$. (We may need to allow $d_1 = 3$ in order to handle singular Toeplitz/Toeplitz-like linear systems; we reduce their solution to computing the rank of the coefficient matrix and to the subsequent solution of a nonsingular Toeplitz-like system.) Like Euclidean algorithm, however, such solution algorithms proceed by recursively decreasing the input size by 1 or by a fixed constant. Therefore, order $N$ parallel steps are required in their implementation. The known alternative algorithms yield NC or randomized NC (RNC) solutions of all the cited computational problems [BGH82], [G84], [BP94], that is, their solution by using $t(n)=O((\log n)^c)$ time and $p(n)=O(n^d)$ arithmetic processors, for 2 constants $c$ and $d$,  under the customary EREW PRAM arithmetic model of parallel computing [KR90], [J92] (cf. our Remark 1.1). (Alternatively, we may define the NC and RNC solutions as the families of arithmetic or aritmetic-Boolean circuits for the above problems having depths $O((\log n)^2)$ and sizes $O(n^d)$ for 2 fixed constants $c$ and $d$ [G86].)Indeed, NC/RNC solution of a linear system of $n$ equations can be computed over any field of constants [Cs74], [Be84], [Ch85], [KP91], [KP92]. There have arisen, however, the important problem of processor efficiency, that is, of having the ratios $p(n)/T(n)$ or even $p(n)/t(n)$ at the level $O((\log n)^{c_1})$ for a constant $c_1$, where $T(n)$ and $t(n)$ denote the record upper and lower bounds on the sequential time of the solution, respectively. Indeed, on the one hand, we have already cited the bound,  $T(n)=O(n (\log n)^3)$, and clearly, $t(n) \geq n$. On the other hand, $p(n)$ has order $n^d$ for $d>3$ in [Cs76], [Be84], [Ch85] and for $d>2$ in [KP91], [KP92], for solving general linear systems in NC/RNC, and $p(n)=O(n^d)$ for $d\geq 2$ for the known NC/RNC Toeplitz/Toeplitz-like solvers over any field of constants [P92], [KP94], [P96]. To yield NC/RNC processor efficiency, we must decrease $d$ to the optimal level 1, and our goal in this paper is to show in detail how this can be done, under certain assumptions on the model of computing. {\bf 1.3. The model of computing.}  Our main assumption is that the input consists of integers (for the reduction from a real input, one may use binary or decimal chopping and then scaling) and that rounding-off a rational number to the closest integer, as well as arithmetic operations, are allowed as unit cost operations. The bit-precision of these computations will also be bounded at the optimal level of the output precision, so that we achieve low Boolean cost. Stating our estimates for the computational cost, we will let $O_A(t,p)$ and $O_B(t,p)$  denote the simultaneous bounds $O(t)$ on the parallel time  and $O(p)$ on the number of arithmetic or Boolean processors, respectively. We will routinely decrease the processor bounds slightly, by exploiting the B-{\it principle} of parallel computing, which is a variant of Brent's principle and according to which $O( s)$ time-steps of a single  arithmetic or Boolean processor may simulate 1 time-step of $s$ arithmetic or, respectively, Boolean processors [KR90], [PP95]. Note that according to the B-principle, the bound $O_A(t, kp)$ implies the bound $O_A(tk, p)$, and similarly $O_B(t, kp)$ implies the bound $O_B(tk, p)$, for a positive parameter $k$. [For $p=1$, we arrive at sequential time bounds $O_A(tk,1)$ and $O_B(tk,1)$.] According to the well-known technique based on the B-principle, one may slow down the computations at the stages requiring too many processors, which in many cases only increases the time bound by a constant factor but more substantially decreases the processor bound ([Q94], pages 44--46, [BP94], pages 297--298). (The more advanced techniques of [PP92], [PP95] sometimes enable an additional increase of processor efficiency but are not needed in our case.){\bf Remark 1.1.} Our algorithms for Toeplitz/Toeplitz-like computations essentially reduce to computing the convolutions (via FFT, assuming that the $2^h$-th roots of 1 are available)  and the inner products of pairs of vectors, and these our basic operations (for vectors of a dimension $n$) can be performed at the cost $O_A(\log n, n)$ and $O_A(\log n, n/\log n)$, respectively, under both EREW PRAM model and more realistic models such as hypercube, butterfly and shuffle-exchange processor arrays [Le92], [Q94]. Thus, it is possible to extend our results to the latter models.{\bf 1.4. Our main results.}Our algorithms of this paper extend our previous work on parallel computations with general matrices [P85], [P87], [P93a], [BP94] (cf. also [PR91], [P92a], [P93b], [PR93]) by means of incorporation of some techniques developed in [P90], [P92], [P92b], [P93], [P93a] for computations with Toeplitz and Toeplitz-like matrices. We will emulate the historic line, by first treating the case of general matrices and then improving the algorithm in the Toeplitz/Toeplitz-like case. We will start with recalling the record parallel complexity bounds of [P85], [P87] for computations with a general $n \times n$ input matrix; we will give  their alternative derivation. Stating these bounds as Theorem 1.1 below, we will use $\omega$ satisfying $2\leq \omega < 2.376$ and such   that a pair of $n \times n$ matrices can be multiplied at the arithmetic cost $O_A(\log n, n^{\omega})$.{\bf Theorem 1.1.} {\it Let $A$ denote a $k \times h$ matrix filled with integers that range from $-2^a$ to $2^a$ for some $a \geq 1$ and for $k +h = O(n)$. Then, with an error probability of at most $n^{-c}$, for a fixed positive constant $c$, one may compute $r$, the rank of $A$, at a randomized computational cost bounded by $O_A((\log n)\log (n a), n^{\omega})$ and by \\ $O_B((\log n) (\log (na)) \log (n\log a) \log \log (n\log a), n^{\omega+1} \log (na))$. If $A$ is nonsingular, then one may compute the inverse and the determinant of $A$ and       solve a linear system $A\vec{x} = \vec{f}$ at a randomized computational cost bounded by $O_A((\log n) \log (na), n^{\omega})$ and by    $O_B((\log n)(\log (n  a))^2 \log \log(n a), n^{\omega+1}\log (\frac{  a}{\log (n a)} +1)\log n)$. The latter cost bounds include the cost of testing correctness of the output.  If $A$ is an $n \times n$ singular matrix, the latter bounds also apply to the computation of an $n \times (n-r)$ matrix formed by $n-r$ basis vectors for a null-space of $A$, and a solution $\vec{x}$ to a linear system $A\vec{f}=\vec{x}$ provided that this system is consitent. The same cost bounds apply to testing correctness of the computed value of rank $A$ as well as of all other output values. In these computations, $2n-1$ random parameters are used for computing $\det A$, and a single random parameter is used for all other tasks. The above complexity estimates do not include the cost of generation of these random parameters.}In the case of an $n \times n$ Toeplitz or Toeplitz-like input matrix, an extension of our approach yields much smaller upper bounds on the complexity of the same computations:{\bf Theorem 1.2.} {\it Let under the assumptions of Theorem~1.1, the $k \times n$ input matrix $A$ be Toeplitz matrix or a Toeplitz-like matrix. Then, all the processor complexity estimates of Theorem~1.1 can be decreased  roughly by factor $n^{\omega -1}/ \log n$, preserving the  time bounds, so as to yield the randomized parallel complexity bounds $O_A((\log n)\log (n a), n \log n)$, \\$O_B((\log n)(\log (na)) \log(n\log a) \log \log (n\log a), n^2 (\log n) \log (na))$,  and   \\    $$O_B((\log n)(\log (n a))^2 \log \log(n a), (\frac{a }{\log ( n a)} +1)(n \log n)^2) \ ,$$ respectively. Here, the inverse of the input matrix and the basis matrix for the null-space of $A$ are assumed to be output in the form of their displacement generators (see Definition 13.2 in our section 13).} {\bf 1.5. Extensions.}We have already cited possible reduction of the computation of each of the {\it gcd}, the {\it lcm}, and the {\it resultant} of 2 polynomials, as well as  {\it Pad\'{e} approximation} of a formal power series or of a polynomial,  to computation of the rank of a Toeplitz/Toeplitz-like matrix and solving a nonsingular Toeplitz/Toeplitz-like linear system of equations (see more details in [BP94], [P96] or [P96,a]). The integrality of the input can be preserved in this reduction, and the input size may grow by a factor of  at most 2. Therefore, {\it the computational complexity estimates of Theorem 1.2 are immediately extended to the listed problems} of the gcd, lcm, resultant and Pad\'{e} computations(assuming the same restrictions on the size and integrality of the input), as well as to various computational problems reducible to the latter ones. Furthermore, we refer the reader to [P90], [H94], and [GKO,a] on the general techniques that immediately enable the extension of our results of Theorem 1.2 to computations with Cauchy-like (generalized Hilbert) and Vandermonde-like input matrices, to [BP93], [BP94] on the extension to the case of matrices represented as the sums of Hankel-like and Toeplitz-like matrices, and to [BGY80], [BP94], [P,b],  and other references cited in the beginning of this paper on various applications of the computations with Toeplitz-like and other structured matrices. Among possible extensions of Theorem 1.1, consider the case where the integer matrix $A$ is symmetric positive definite, sparse, and associated with an $s(n)$-separatable graph given with its $s(n)$-separator family (cf. [LRT79], [P93b], [PR93]).  If such a matirx $A$ is well-conditioned (even if its entries are not integers but any real numbers), then, at the arithmetic cost $O((\log n)^3, (s(n))^{\omega}/\log n)$, the parallel algorithm of [P93b], [PR93] numerically computes both recursive factorization of such a matrix and its determinant, as well as a solution $\vec{x} = A^{-1} \vec{f}$ to a linear system $A\vec{x}=\vec{f}$ (if $\det A \neq 0$). Numerical approximation is involved in this algorithm at the auxiliary stages of matrix  inversions, where a parallel algorithm of [PR89] is applied. If $A$ is filled with integers, then this stage can be performed exactly, by using the algorithms of [P85], [P87]. Then, the exact recursive factorization of $A$, $\det A$, and $A^{-1}\vec{f}$ can be computed at the arithmetic cost $O_A((\log n )^3, (s(n))^{\omega} +n)$.  By employing the algorithm of this paper for general matrix inversion, one may improve the latter bounds a little, so as to yield $O_A((\log n )^2, (s(n))^{\omega} +n)$. The results of Theorems 1.1 and 1.2 can be further extended to various other matrix computations, by using the known reduction techniques of [BP94], [P96], [P96b]. For demonstration, consider the computation of the characteristic polynomial $c_A(x)=\det (xI-A)$ of the above sparse $n \times n$ matrix $A$. Such a polynomial has degree $n$. We may first concurrently compute $c_A(x)$ at $n+1$ distinct points $x_0,\ldots, x_n$ and then interpolate so as to obtain the coefficients of $c_A(x)$. If the chosen values of $x_i$ are much larger than $\|A\|$, then the matrices $x_iI - A$ are positive definite, and we may compute $c_A(x_i)$ for $i=0, 1, \ldots, n$, at the overall computational cost $O_A((\log n )^2, (s(n))^{\omega} +n)n)$. These bounds dominate the cost of the subsequent interpolation producing the polynomial $\det (xI -A)$. {\bf 1.6. Some technical and historical comments.} Historically, as we have mentioned already, Theorem~1.1 was proved in [P85], [P87]. Furthermore, the techniques of [P85] and [P87] are the springboard for our present construction, and our reliance on these techniques is the reason for choosing the model of computations with integer input and with rounding-off allowed at unit cost. The main techniques of [P85] and [P87], extensively used in the present paper (as well as in [BP94], p. 357, and [R95]), include combining the algebraic and numerical versions of Newton's iteration for matrix inversion, $p$-adic (Newton-Hensel's) lifting, and the variable diagonal method introduced in [P85], [P87] in order to relate the above techniques to each other  (see more comments in our Appendix D). A slightly weaker version of Theorem~1.2 (with time-bound larger by factor $n$) has been deduced in [BP94], p. 357, by means of combining the approach of [P85], [P87] with the well-known techniques for recursive block factorization of a matrix [St69], [BA80], [Morf80], [PR91], [P93b], [PR93] and for Toeplitz-like computations, including, in particular, more recent techniques of [P92], [P93], and [P93a], in particular the  {\it techniques of the SVD truncation of a displacement generator}. In [R95] the bounds $O((\log n)^2, (\log n)^{\omega}n)$, for $2\leq \omega <2.376$ (improving the arithmetictime-bound of [BP94] by factor $\log n$, though increasing by  factor $(\log n)^{\omega}$, $2\leq \omega <2.376$, the arithmetic processor bound of [BP94]) and the bound $O(n \log n)$ on the bit-precision of the computations have been claimed under the additional restriction that all the entries of the input matrix were integers  ranging  between $-n^{\beta}$ and $n^{\beta}$ for a constant $\beta$ independent of $n$. The algorithm of [R95] is based on the idea of incorporation of the pipelining techniques of stream contraction of [PR91] into the constructions that follow the lines of either [BP94] or [P85], [P87], [BA80] and [Morf80]. In fact, such a simple and nice idea is the only novelty of a long paper [R95]. As we can see from the present paper, this idea has a nontrivial implementation that gives us the desired minor improvement of the time bound of [BP94]. The algorithm of [R95], however, was presented with at least 3 major technical errors (see Remarks~8.1, 11.1, 14.1, and 14.2 in our sections 8, 11, and 14). Due to one of these errors (see our Remark~14.1), the algorithm of [R95] falls apart. Even if this serious error has been fixed, then, due to the 2 other errors,    the algorithm of [R95] would have only supported  upper bounds on the numbers of arithmetic and Boolean processors exceeding the bounds claimed in [R95] by factors at least $n^{\omega -1 }/ \log n$, $2 \leq \omega < 2.376$. To fix all these 3 errors, one needs to involve some nontrivial techniques, which are much more sophisticated and advanced (both technically and conceptually) than the cited simple nice idea of [R95].{\bf Remark 1.1.} Apart from the 3 cited major errors, the presentation in [R95] is still technically careless. For instance, the main algorithm of [R95] uses randomization,  but this is not mentioned in the statements of the resulting complexity estimates. In another instance, the upper estimates of Lemma 2.4 on page 753 of [R95] are fundamental for the main algorithm of [R95], but their proofs are given neither in [R95] nor in the bibliography cited in conjunction  with this lemma. In fact, part of this bibliography  is irrelevant to Lemma 2.4. Moreover, the upper estimates stated in the part (b) of the lemma  are much greater than the known ones, whereas the upper estimates stated in the part (c) are " better" (that is, less) than the record estimates presently available in the literature. Furthermore, besides the 3 cited major errors, the paper [R95] contains several flaws in the proofs and less desastrous  but still serious errors, as well as rediscoveries of the results and techniques of [P85], [P87], and [BP94] and irrelevant and/or invalid citations. Incidently, [R95] cites the paper [R93], as being related to [R95], and [R93] itself also has major drawbacks. In particular, [R93] claimed a dramatic improvement (by factor $n$) of the known complexity results, but the comparison of [R93] with [BT90] and [BP91] shows immediately that the claim should be narrowed from factor $n$ to at best factor $\log n$. Moreover, unlike the algorithms of [BT90] and [BP91] and like the algorithm of [R95], the algorithm of [R93] is incorrect or at least does not support the claimed complexity bounds. A disturbing feature of [R95] is incorporation of almost all the major techniques of [P85], [P87] (including introduction of a variable diagonal, rapid approximation via Newton's iteration, rounding-off approximations to the closest integer, reduction modulo a prime, Hensel's lifting, and the final recovery of the output from its value modulo a prime power) in such a way that the novice reader should attribute all these techniques to [R95]. Such a reader may also exaggerate the value of the paper [R95] by believing, for instance, in Reif's claim about the extension of the results of [R95] to BCH decoding. This claim is deceptive because the application to coding is over finite fileds, where Reif's claimed results would have being highly ineffective even if he could have supported all his technical claims. In our present paper, we give a detailed proof of Theorems~1.1and 1.2, by revisiting the approach of [BP94], by supplying several omitted details, and by incorporating the stream contraction techniques of [PR91]. The resulting algorithm (supporting Theorem~1.1) elaborates the block Gauss-Jordan decomposition of a matrix filled with integers, so as to produce the output values as rationals with bounded denominators; then, such rationals are made integers via scaling. The techniques that ensure the bounds on the denominators are novel and may be of independent interest. They extend (from the case of Gaussian elimination to the recursive block Gauss-Jordan case) the celebrated techniques of symbolic solution of linear systems [E67], [B68]. To improve parallel computation of the recursive block Gauss-Jordan decomposition, we perform the auxiliary steps of matrix inversion approximately, by means of Newton's iteration,  for which we develop a particularly effective parallel implementation in the Toeplitz and Toeplitz-like cases. The resulting approximation is then extended to the exact rational decomposition, which is computed, in particular, by means of rounding-off the computed approximations to the closest integers. The initial approximations for Newton's iteration,  applied in order to perform all the auxiliary matrix inversions, is ensured by means of varying the diagonal of the input matrix, that is, shifting to strongly diagonally dominant input matrix $A+psI$, for a fixed pair of a prime $p$ and an integer $s$ (cf. [P85], [P87]). The decomposition of $A+psI$ gives us the decomposition modulo $p$ of $A$, which we lift (via Newton-Hensel's lifting) so as to arrive at the decomposition modulo $p^{2^g}$ of $A$, for  $g=\lceil \log ((n \log \|A\|)/\log p)\rceil +1$. By using the latter decomposition, we compute $(\det~A) \bmod p^{2^g}$ and $({\rm adj}~A) \bmod p^{2^g}$. Then, since $g$ is large enough and since $\det~A$ and the entries of ${\rm adj}~A$ are integers, we easily recover them. Then we compute $A^{-1}={\rm adj}~A/\det~A$ if $\det~A \neq 0$. Such a construction and the various algebraic and numerical techniques involved may be of independent theoretical interest. {\bf 1.7. Organization of the paper.}We will organize our paper as follows. After some preliminaries in section~2, we will introduce recursive decomposition (RD) and extended RD (ERD) of a matrix in section~3. We define integral RD (IRD) in section~4, show a technique for transition from RD to IRD in section~5, and estimate parallel complexity of computing the IRD of a diagonally dominant matrix in section~6. We recall an algorithm for approximate matrix inversion via Newton's iteration in section~7, apply it in order to approximate the RD of a matrix, with extension to computing the IRDof an integer matrix, in section~8, and estimate parallel complexity of these computations in section~9. We apply pipelining (stream contraction) to achieve acceleration by factor $\log n$ in section~10, extend the results of sections 8 and 10 to computing the ERD modulo a fixed prime in section~11, and use $p$-adic lifting to recover from this ERD the inverse, the determinant, the rank, and the null-space of an integer matrix (thus proving Theorem~1.1) in section~12. In section~13, we recall some known definitions and properties for computations with Toeplitz and Toeplitz-like matrices, and in section~14, we apply these properties in order to improve the results of section~12 in the Toeplitz and Toeplitz-like cases. Section~15 is left for discussion. Some details of our proofs are shown in the parts A, B, and C of the Appendix, whose part D contains a brief recollection of some of the main techniques of [P85], [P87], incorporated in the algorithms of the present paper. \vspace{.23in}%\noindent\begin{center}{\bf 2. Some Definitions and Auxiliary Results}\end{center}{\bf Definition 2.1.} $I$ and $O$ denote the identity and  null matrices, of appropriate sizes. $W^T$ is the transpose of a matrix or vector $W$.  diag $(w_i )_{i=0}^{n-1}=$ diag $(w_0 ,\ldots ,w_{n-1})$ is the diagonal matrix whose diagonal isfilled with $w_0 ,\ldots ,w_{n-1}$; $D(W)=$ diag $(W)=$ diag$(w_{i,i})_{i=0}^{n-1}$ for a matrix $W=(w_{i,j})$.rank $W$ is the rank of $W$. det $W$ is the determinant of a square matrix $W$. adj $W$ is the adjoint (adjugate) matrix of $W$, equal to $W^{-1}$ det $W$ for a nonsingular matrix $W$.{\bf Definition 2.2} [GL89].$\parallel {\bf v}\parallel= \parallel {\bf v}\parallel_1 =\displaystyle{\sum_i |v_i |}$,$\parallel {\bf v}\parallel_2= \displaystyle{\left(\sum_i v_i^2 \right)^{1/2} }$,for a real vector ${\bf v} =(v_i )$.$\parallel W\parallel_g= \max_{\parallel v\parallel_g =1} \parallel Wv\parallel_g$,$g=1,2$; $\parallel W\parallel =\parallel W\parallel_1$,for a matrix $W$.{\bf Proposition 2.1} [GL89]. $\parallel W\parallel =\parallel W\parallel_1 = \max_j \displaystyle{ \sum_i |w_{i,j} |}$,{\it for a matrix $W=(w_{i,j})$, and if $W$ is a $k\times k$ matrix, then}$$\parallel W\parallel /k^{1/2} \leq\parallel W\parallel_2 \leq \parallelW\parallel k^{1/2} \ .$${\bf Definition 2.3.}$d(W) =\parallel WD^{-1}(W)-I\parallel$. A matrix $W$ is{\it column-diagonally dominant (c.-d.d.)} if $d(W) <1$.{\bf Definition 2.4.} For a $k\times k$ matrix $W$, let$W^{(q)}$ denote its $q\times q$ {\it leading principal submatrix(l.p.s.),} formed by the intersection of the first $q$ rows and thefirst $q$ columns of $V$, $q=1,2,\ldots ,k$. If $B=W^{(q)}$ is the$q\times q$ nonsingular l.p.s. of $W$ and if$$W =\left( \begin{array}{cc}B & C \\ E & G \end{array}\right) , \eqno (2.1)$$then the matrix$$S = S(W,B) = G-EB^{-1}C \eqno (2.2)$$is called the {\it Schur complement} of $B$ in $W$. By applying blockGauss-Jordan elimination to the $2\times 2$ block matrix $W$ of (2.1),we obtain{\bf Proposition 2.2.} {\it Under the assumptions of Definition 2.4,we have$$W =\left( \begin{array}{cc}I & O \\ EB^{-1} & I \end{array}\right) \ \left( \begin{array}{cc}B & O \\ O & S \end{array}\right) \ \left( \begin{array}{cc}I & B^{-1}C \\ O & I \end{array}\right) \ , \eqno (2.3)$$and if the matrix $W$ is nonsingular, then also $S$ is nonsingular, and  }$$W^{-1} =\left( \begin{array}{cc}I & -B^{-1}C \\ O & I \end{array}\right) \ \left( \begin{array}{cc}B^{-1} & O \\ O & S^{-1} \end{array}\right) \ \left( \begin{array}{cc}I & O \\ -EB^{-1} & I \end{array}\right) \ . \eqno (2.4)$$We also have the following easily verified results:{\bf Proposition 2.3} (cf. [BP94]). {\it For $k>q$, the Schur complement of$B=W^{(q)}$ in a $k\times k$ matrix $W$ is the $(k-q)\times (k-q)$matrix obtained from $W$ in $q$ steps of Gaussian elimination(without pivoting) provided that these $q$ elimination steps can becarried out (with no division by 0). The latter assumption holds, inparticular, if $d(W)<1$.}{\bf Proposition 2.4}, transitivity of Schur's complementation. {\it If $W^{(r)}$ and $W^{(q)}$ are nonsingular matrices, then  $S(W,W^{(q)})=S(S(W,W^{(r)}), S(W,W^{(r)})^{(r-q)})$.} {\bf Corollary 2.5.} {\it If $d(W)<1$ for a matrix $W$ of (2.1), then$B$, $S$, and $W$ are nonsingular matrices,$d(B)\leq d(W) <1$, $d(S) \leq d(W)<1$, and(2.3)--(2.4) hold.}{\bf Corollary 2.6.} {\it If a matrix $W$ satisfies (2.1)--(2.3) andis filled with integers, then so is the matrix $S$ {\rm det}$B$.}{\bf Proposition 2.7.} {\it The matrix equation (2.3) implies that{\rm det} $W=$ ({\rm det} B) {\rm det} $S$, and furthermore, that $$W^{(t)} =\left( \begin{array}{cc}I & 0 \\ EB^{-1} & I \end{array}\right)^{(t)} \ \left( \begin{array}{cc}B & O \\ O & S \end{array}\right)^{(t)} \ \left( \begin{array}{cc}I & B^{-1}C \\ O & I \end{array}\right)^{(t)},\ \ \ {\rm for} \ t=1,\ldots ,k  \ , $$and det $W^{(t)} =$(det $B$) det $S^{(t-q)}$, for $t=q+1,\ldots k$, provided that $W$ is a $k \times k$ matrix.}We will also use the following known fact (cf. [GL89] or [BP94]:{\bf Proposition 2.8.} {\it$|{\rm det} \ W|\leq (\parallel W\parallel_g)^k$, and furthermore, $|v| \leq (\parallel W \parallel_g)^{k-1}$ for every entry $v$ of adj $W$, for $g=1,2$ and for a $k\times k$ matrix $W$.}{\bf Definition 2.5.} A matrix $W$ is {\it strongly nonsingular} ifall its leading principal submatrices are nonsingular.{\bf Definition 2.6.} A real matrix $M$ is {\it symmetric positivedefinite (s.p.d.)} if it can be represented as the product $A^T A$ fora nonsingular matrix $A$.{\bf Proposition 2.9.} (cf. [GL89], [BP94], exercise 4c, page212). {\itIf a matrix $W$ of (2.1) is strongly nonsingular (respectively, iss.p.d.), then so are its every l.p.s. [including the matrixB of (2.1)] and the Schur complement S of B, defined by (2.2).}{\bf Corollary 2.10.} {\itAny s.p.d. matrix is strongly nonsingular.}{\bf Proposition 2.11.} {\it For a real matrix $A$ and the matrix $AA^T$, we have ${\rm rank}~A = {\rm rank}~(AA^T)= \max \{ k \ \  such \ \  that \ \  \det (AA^T)^{(k)} \neq 0 \}$.} In section 12, we will also use the following definition and simple result: {\bf Definition 2.7} (cf. [GL89] or [BP94]). {\it The null-space ${\bf N}(A)$ of a matrix $A$ is the linear space formed by all vectors $\vec{x}$ satisfying the vector equation $A\vec{x}=\vec{0}$.} {\bf Fact 2.12.} {\it 2 vectors equations, $A\vec{x}=\vec{f}$ and $A\vec{y}=\vec{f}$, together imply that $\vec{x} - \vec{y} \in {\bf N}(A)$ or, equivalently, any solution $\vec{x}$ to a consistent linear system $A\vec{x}=\vec{f}$ can be represented in the form $\vec{x} = \vec{x_0}+\vec{z}$, where $\vec{x_0}$ is a fixed specific solution and $\vec{z} \in {\bf N}(A)$.}\vspace{.23in}%\noindent\begin{center}{\bf 3. Recursive Decomposition (RD)and Extended RD (ERD)of a C.-D.D. Matrix}\end{center}Hereafter, for convenience, let $\log$ stand for $\log_2$, let $n=2^h$ for an integer $h= \log n$,and let $V$ be a given $n\times n$ c.-d.d. matrix. We will define a{\it recursive decomposition} {\it (RD)} of this matrix based on itsrepresentation in the form (2.3)for $q=n/2$; we will first apply (2.3) to$W=V$ and then, recursively, to $W=B$ and $W=S$, and so on, though in fact, we willmostly care about the diagonal blocks $V,V_0 =B$, $V_1 =S$, and so on, whichwe will identify with the nodes of a binary tree, $T$.The node of $T$ associated with a binary string $\alpha$ of length $|\alpha |$ is a $k\times k$ matrix, denoted $V_{\alpha}$, where$k=n/2^{|\alpha |}$. The root of the tree is the $n\times n$ matrix$V=V_{\Lambda}$, associated with the empty string $\Lambda$. For a string $\alpha$, we let$\alpha 0$ and $\alpha 1$ denote the 2 strings obtained by appending 0 and1 to $\alpha$, respectively. Then, we will assume that the matrixequations (2.1)--(2.4) are satisfied for $W=V_{\alpha}$,$B=V_{\alpha 0}$, $S=V_{\alpha 1}$, and for any binary string $\alpha$of length $|\alpha | <h$. The resulting RD of the matrix $V$ continuesup to the level $h$, where it reaches its leaves-matrices$V_{\alpha} =(v_{\alpha} )$ of size $1\times 1$, where$|\alpha |=h$, and  where $v_{\alpha}$ denotes the single entry of$V_{\alpha}$.{\bf Proposition 3.1.} {\it All the nodes $V_{\alpha}$ of the tree $T$are c.-d.d. matrices; moreover, $d(V_{\alpha}) \leq d(V) < 1$ forall strings $\alpha$.}{\bf Proof.} Recall that $V$ is an c.-d.d. matrix and recursively applyCorollary 2.5.\hfill $\Box$Let us formalize the computation of the RD of $V$, by extendingthe notation (2.1)--(2.2) to $V_{\alpha}$. We will write$$V_{\alpha} = \left( \begin{array}{cc} B_{\alpha} & C_{\alpha} \\E_{\alpha} & G_{\alpha} \end{array} \right) \ , \eqno (3.1)$$$$V_{\alpha 0} = B_{\alpha} \ , \ \ \V_{\alpha 1} =S_{\alpha} = G_{\alpha} -E_{\alpha} B_{\alpha}^{-1}C_{\alpha} \ , \eqno (3.2)$$for all binary strings $\alpha$ of length less than $h$.Then, computation of the RD of a c.-d.d. matrix $V$ amounts to recursivecomputation of $V_{\alpha 1}$ of (3.2) for all binary strings $\alpha$of length increasing from 0 to $h-1$. The computation produces thematrices $V_{\alpha 0}^{-1} =B_{\alpha}^{-1}$ for all binary strings$\alpha$ of length less than $h$, as auxiliary results. By appending theseinverse matrices $V_{\alpha 0}^{-1}$ to the nodes $V_{\alpha 0}$ of thetree $T$ (and, consequently, to the RD of $V$), we arrive at the{\it extended RD} of $V$, which we will abbreviate as the ERD of$V$. The set of the matrices $V_{\alpha 0}^{-1}$ will be called the{\it extending set of the RD of $V$}.The RD and the ERD of any matrix $V$ can be defined as long as all theinvolved nodes-matrices $V_{\alpha 0}$ for all the binary strings$\alpha$ of length less than $h$ are nonsingular. Recursive applicationof Proposition 2.9 and Corollary 2.10 yields{\bf Proposition 3.2.} {\itThere exists the RD and the ERD of any strongly nonsingular (inparticular, of any s.p.d.) matrix.}{\bf Remark 3.1.} As soon as we have the RD of a c.-d.d. matrix $V$,we immediately obtain the RD of $V^{-1}$ based on recursive applicationof (2.4). By using such an RD, we may compute the solution$\vec{x} = V^{-1} \vec{f}$ of a linear system $V\vec{x}=\vec{f}$, at alower computational cost $O_A((\log n)^2, n^2/(\log n)^2)$. {\bf Remark 3.2,}{\it an overview of the next sections.}We will need to compute the RDs of c.-d.d. integer matrices $V$ of some special class (see section 9), which then will enable us to compute $\det V_{\alpha}$ and $\det V^{(q)}$ for all binary strings $\alpha$ of length at most $h$ and for $q=2, \ldots, n$, at a lower computational cost $O_A(\log n, n/\log n)$(see Propositions 5.1 and 5.2 in section 5). This will be an auxiliary steptowards devising an efficient parallel algorithm for computation of thedeterminant and the inverse of any nonsingular matrix filled withintegers. The algorithm will be shown to be particularly effectivein the case of a Toeplitz or Toeplitz-like input matrix (filled withintegers).\vspace{.23in}%\noindent\begin{center}{\bf 4. Integral Recursive Decomposition (IRD) of a C.-D.D. MatrixFilled with Integers} \end{center}Suppose that the matrix $V$ is filled with integers. Then, for allbinary strings $\alpha$, the matrices $V_{\alpha}$ are filledwith rationals, and there exist integer multipliers $m_{\alpha}$ suchthat $m_{\alpha} V_{\alpha}$ are integer matrices.We will next specify a particular choice of such integer multipliers$m_{\alpha}$.We will use the notation $W^{(q)}$ of Definition 2.4 and the following definition:{\bf Definition 4.1.} $(\alpha)_2$ denotes the binary valuerepresented by a binary string $\alpha$ (of length at most $h$).$\alpha (q)$ denotes the binary string that represents a nonnegative binaryvalue $q$, so that $(\alpha (q))_2 =q$ and $\alpha ((\alpha )_2 )=\alpha$.$H(|\alpha |)$ denotes $2^{h-|\alpha |}=n/2^{|\alpha|}$. $Q(\alpha )$ denotes $(\alpha)_2 H(|\alpha |)$.By applying Propositions 2.3 and/or 2.4, we obtain{\bf Proposition 4.1.} {\it Let a binary string $\alpha$ end withbit 1 and have length at most $h$. Then the matrix $V_{\alpha}$is the Schur complement of$V^{(Q(\alpha ))}$ in $V^{(Q(\alpha) + H(\alpha))}$. }By combining this proposition with Corollary 2.6, we obtain{\bf Corollary 4.2.} {\it If the matrix V is filled with integers, thenso are the matrices $m_{\alpha 1}V_{\alpha 1}$, for $m_{\alpha 1} ={\rm det} \ V^{(Q(\alpha 1))}$ and for all binary strings $\alpha$ oflength less than h}.By replacing the matrix $V_{\alpha 1\gamma}$ bythe pair of the scalar $m_{\alpha 1}= m_{\alpha 1 \gamma}$ and the matrix $m_{\alpha 1} V_{\alpha 1\gamma}$, in the binary tree $T$, for every pair of binary strings$\alpha$ and $\gamma$ such that $|\alpha |+|\gamma |<h$ and $\gamma$consists only of zeros, we arrive at the {\it integral recursivedecomposition} (IRD) of a c.-d.d. integer matrix $V$. By definition, we will also include the set $\{ \det A^{(k)}, \  k=1,\ldots, n\}$ into the IRD of $A$. Clearly, having the IRD of $V$ available, we may immediatelty compute the RD of $V$. In the next section, we will show how to obtain the IRD from the RD.{\bf Remark 4.1.} Due to Proposition 2.8, we have$$|m_{\beta \gamma}| =|{\rm det} \ V^{(Q (\beta ))}|\leq \parallel V^{(Q(\beta ))}\parallel^{Q(\beta )} \ , $$ if  the binary string  $\beta$ ends  with  bit 1 and if $\gamma$ is  a string  of  zeros, and therefore,$$|m_{\alpha}| \leq \parallel V\parallel^n \ , \ \ \ {\rm for \ all \ binary \  strings \   } \alpha, \ |\alpha| \leq h \ . \eqno (4.1)$$\vspace{.23in}%\noindent\begin{center}{\bf 5. Transition from RD to IRD}\end{center}For a c.-d.d. integer matrix $V$ given together with itsRD, we will seek its IRD. We recall that $v_{\alpha}$ for$|\alpha |=h$ denotes the single entry of the $1\times 1$ leaf-matrix$V_{\alpha}$ of the tree $T$ of the RD. Recursive application ofProposition 2.7 yields the two following results:{\bf Proposition 5.1.} {\it For every binary string $\alpha$ oflength at most h, we have$${\rm det} \ V_{\alpha} = \prod_{\beta } v_{\beta} \ ,$$where $\displaystyle{ \prod_{\beta } }$ denotes theproduct in all binary strings $\beta$ of length h that have $\alpha$ as their prefix, that is, the associated nodes $V_{\beta}$ are both leaves of the tree $T$ and descendants of the node $V_{\alpha}$ in the tree T}.{\bf Proposition 5.2.} {\it $\det~V^{(q)} =  \displaystyle{\prod_{(\alpha)_2 < q }}$$v_{\alpha}$, where $\displaystyle{\prod_{(\alpha)_2 <q} }$ denotes the product in all binary strings$\alpha$of length $h$ for which $(\alpha )_2 <q$.}{\bf Definition 5.1.} $\det V^{(q)}$ and $\det V_{\alpha}$, for all $q$ and $\alpha$, will be called  {\em associated to the RD} of $V$. From Propositions 4.1, 5.1, 5.2, and Corollary 2.1, we obtain the following result, by applying the well-known parallel prefix algorithm [EG88], [KR90]:{\bf Corollary 5.3.} {\it Given the RD of a c.-d.d. $n\times n$matrix $V$ filled with integers, one may compute the IRD of V and all the associated determinants of the RD of $V$, at thecost $O_A(\log n,n/\log n)$.}\vspace{.23in}%\noindent\begin{center} {\bf 6. Complexity of Computing the Integral Recursive Factorization(IRD) of a C.-D.D. Matrix Filled with Integers} \end{center}Let a c.-d.d. $n\times n$ matrix $V$ be filled with integersand be given together with the integer values of $|$det $V^{(q)}|$ for$q=2,3,\ldots ,n$. Then, due to the matrix equations (2.1)--(2.3)and to the results of the previous section, the computation of the IRDof $V$ can be reduced essentially to a sequence of multiplications,inversions, and subtractions of integer matrices. We recall the knowncomplexity estimates:$$O_A(t_S (k), p_S (k))=O_A(1,k^2 ) \eqno (6.1)$$for $k\times k$ matrix subtraction (obvious),$$O_A(t_M (k),p_M (k))=O_A(\log k,k^{\omega} ), \ \ \\ 2 \leq \omega < 2.376 \ , \eqno (6.2)$$for $k\times k$ matrix multiplicaton [BP94], and$$O_A(t_I (k),p_I (k))=O_A((\log k)\log (k\log \parallel W\parallel ),k^{\omega})\ , \ \ \ 2 \leq \omega < 2.376 , \eqno (6.3)$$for the inversion of a $k\times k$ nonsingular matrix $W$ filled withintegers [P85], [P87] (provided thatusing a single random integer parameter is allowed forperforming the inversion).Let us write$$v=\max_{\alpha} \{\parallel V_{\alpha} \parallel\} \ , \eqno (6.4)$$recall Remark 4.1, and conclude that we may substitute$$\parallel W_{\alpha} \parallel \leq v\parallel V\parallel^n \eqno (6.5)$$in (6.3),as long as $W_{\alpha} =m_{\alpha} V_{\alpha}$ for a binary string $\alpha$,$|\alpha |\leq h$.Now, we will extend the above relations (6.1)-(6.5) so asto estimate the overall complexity of computing the IRD of $V$.The matrix equations (2.1) --(2.3) imply that$$t_{\rm IRD} (n)\leq t_I (n/2)+t_{\rm IRD} (n/2) +2t_M (n/2) + 1 ,$$$$p_{\rm IRD}(n)\leq \max \{p_I (n/2), 2p_{\rm IRD} (n/2), p_M (n/2), n^2 \} \ .$$By recursively applying the latter bounds, together with the bounds (6.1)-(6.5), and by using the B-principle, which we have recalled in section 1, we arrive at the following result:{\bf Proposition 6.1.} {\it The IRD of a c.-d.d. $n\times n$matrix V filled with integers can be computed at the parallel arithmetic cost bounded by $O_A((\log n)^2 \log ((\log v)+ n \log \parallel V\parallel )),n^{\omega}/\log n)$,for $\omega $ of (6.2),(6.3) and for v of (6.4).}\vspace{.23in}%\noindent\begin{center}{\bf 7. Approximate Matrix Inversion via Newton's Iteration}\end{center}In sections 9--12, we will modify and refine the stage of the auxiliary matrix inversions in the algorithm supporting proposition 6.1, which will enable us to improve its estimates. In section 13 and 14, we will further adjust the construction to the case of Toeplitz/Toeplitz-like input. In this section we will recall approximate matrix inversion via Newton's iteration. Then, in the next sections, we will modify the construction of section 6, so as to include this iteration and, later on, to recover the exact inverses, by means of rounding-off rationals to the closest integers.{\bf Algorithm 7.1}, {\it Newton's iteration for approximate matrixinversion.}{\bf Input:} a nonsingular $k\times k$ matrix $B$, two positive scalars$b$ and $c$, $b>c$, and a matrix $X_0$ (of a rough initialapproximation to $-B^{-1}$) such that$$c=-\log\parallel BX_0 +I\parallel \ . \eqno (7.1)$${\bf Output:} a matrix $X$ such that$$\parallel BX-I\parallel \leq 2^{-b} \eqno (7.2)$$and, consequently,$$\parallel B^{-1}-X\parallel\leq 2^{-b}\parallel B^{-1}\parallel \ .\eqno (7.3)$$\begin{enumerate}\item Compute$$g=\lceil \log (b/c)\rceil \ . \eqno (7.4)$$\item Recursively compute the matrices$$X_i =X_{i-1}(2I+BX_{i-1}) \ , \ \ \ \ i=1,\ldots ,g \ . \eqno (7.5)$$\item Output the matrix $X=-X_g$.\end{enumerate}To prove {\it correctness} of Algorithm 7.1, deduce from (7.5) that$$I+BX_i=(I+BX_{i-1})^2 \ , \ \ \ \ i=1,2,\ldots ,g \ ,$$and consequently,$$I+BX_i =(I+BX_0 )^{2^i} \ ,$$$$\parallel I+BX_i \parallel\leq\parallel I+BX_0 \parallel^{2^i} \ ,\eqno (7.6)$$for $i=1,2,\ldots ,g$. In particular, for $i=g$, we have$$\parallel I+BX_g \parallel\leq\parallel I+BX_0 \parallel^{2^g}=2^{-c2^g} \leq 2^{-b} \ ,$$due to (7.1) and (7.4). Thisgives us (7.2) and (7.3), for $X=-X_g$.To estimate the overall {\it computational cost}  of performing Algorithm 7.1,observe that the $i$-th step (7.5) amounts to 2 matrixmultiplications and to adding, at the cost $O_A(1,k)$, the matrix$2I$ to the matrix $BX_{i-1}$. Summarizing, we obtain the followingresult:{\bf Proposition 7.1.} {\it For g of (7.4), g steps (7.5) of Newton'siteration, performed at the overall cost $O_A(gt_M (k),p_M (k))=O_A(g\log k,k^{\omega} )$, for $\omega$ of (6.2),   $2\leq \omega <2.376$, suffice in order to compute a matrix$X=-X_g$ satisfying (7.2) and (7.3).}{\bf Remark 7.1.} If $B$ is a c.-d.d. matrix, so that $d(B) <1$, then(7.1) holds for $X_0 =-D^{-1} (B)$ and $c=-\log d(B)$.\vspace{.23in}%\noindent\begin{center}{\bf 8. Approximate RD of an Integer Matrix and Its Extension to theExact Evaluation of the IRD}\end{center}Recall the algorithm of sections 3--6 for computing the ERD, IRD, and the associated determinants of a c.-d.d.integer matrix $V$ but change its blocks performing matrix inversions soas to use Algorithm 7.1 for some fixed sufficiently large $b$ and for$X_0$ and $c$ defined according to Remark 7.1. Then, the matrices$V_{\alpha}$ and the values of det $V^{(q)}$ and det $V_{\alpha}$ for all$q$ and $\alpha$ are computed approximately, even though we stillperform all arithmetic operationsover the rationals, with infinite precision andno errors.Let us write$\tilde{V}_{\alpha}$, $\tilde{\det}\ V^{(q)}$, and  $\tilde{\det}\  V_{\alpha}$ for the computed approximations to $V_{\alpha}$, $\det \ V^{(q)}$, and $\det \ V_{\alpha}$, respectively. Then, we extend the relations(3.1), (3.2), (7.1), (7.4), and (7.5) by writing$$\tilde{V}_{\alpha} = \left( \begin{array}{cc}\tilde{B}_{\alpha} & \tilde{C}_{\alpha } \\\tilde{E}_{\alpha} & \tilde{G}_{\alpha } \end{array}\right) \ ,\eqno (8.1)$$$$\tilde{V}_{\alpha 0} = \tilde{B}_{\alpha} \ ,\ \ \ \ \ \tilde{V}_{\alpha 1} =\tilde{G}_{\alpha}-\tilde{E}_{\alpha} X_{\alpha} \tilde{C}_{\alpha} \ , \eqno (8.2)$$for all binary strings $\alpha$ of length less than $h$, where $$X_{\alpha ,0} =-D^{-1}(\tilde{B}_{\alpha} ) \ , \eqno (8.3)$$$$c(\alpha )=-\log\parallel\tilde{B}_{\alpha}X_{\alpha ,0}+I\parallel \ , \eqno (8.4)$$$$g(\alpha )=\log (b/c(\alpha )) \ , \eqno (8.5)$$$$X_{\alpha ,i}=X_{\alpha ,i-1} (2I+\tilde{B}_{\alpha} X_{\alpha ,i-1} ) \ , \ \ \ \ i=1,\ldots ,g(\alpha ) \ , \eqno (8.6)$$$$X_{\alpha} =-X_{\alpha ,g(\alpha )} \ .  \eqno (8.7)$${\bf Algorithm 8.1,} {\it approximating the RD of a c.-d.d. matrix}.{\bf Input:} a positive integer $h$, $n=2^h$, a positive $b$, and an $n\times n$ c.-d.d. matrix $V=\tilde{V_{\Lambda}}$.{\bf Output:} a set of matrices $\tilde{V}_{\alpha}$ for all binarystrings $\alpha$ of length at most $h$, satisfying (8.1)--(8.7), and the values $ \tilde{\det}~V_{\alpha}$ for all $\alpha$ and  $\tilde{\det}~V^{(q)}, \ q=1, \ldots, n$, defined according Proposition~5.2 with det and $v$ replaced by $\tilde{\det}$ and $\tilde{v}$, respectively.{\bf Computations:} Recursively apply the equations (8.1)--(8.7),starting with $V=\tilde{V_{\Lambda}}$, for the empty string $\Lambda$, in order to compute$\tilde{V}_{\alpha 0}$ and $\tilde{V}_{\alpha 1}$ for all binarystrings $\alpha$ of length less than $h$. (For binary strings $\alpha$of length $h-1$, the matrices $\tilde{V}_{\alpha 0}$ and $\tilde{V}_{\alpha 1}$ have size $1\times 1$, so$\tilde{V}_{\alpha 0}$ is inverted immediately, and thenthe matrix $\tilde{V}_{\alpha 1}$ is computed based on (2.1) and (2.2) for $W=\tilde{V_{\alpha}}$.) Finally, compute ${\rm \tilde{\det}}~V_{\alpha}$ and  ${\rm \tilde{\det}}~V^{(q)}$ for all $\alpha$ and $q$ by applying Propositions 5.1 and 5.2, under the above modification of the notation.{\it Correctness} of the algorithm is immediately verified. {\it The computational cost} of performing it is bounded by $O_A((\log n)^2 g, n^{\omega})$, for $g = \max_{|\alpha| \leq h} g(\alpha)$ and for $\omega$ of (6.2). [Repeat the proof of Proposition~6.1 but rely on Proposition~7.1 instead of (6.3).]Next, for a c.-d.d. matrix $V$ filled with integers, we willapply Algorithm 8.1 for a sufficiently large $b$, and then, we  will apply the techniques of {\it integer rounding-off} (compare [P85], [P87], and [BP94], page 252) in order to extend the resulting approximate RD of$V$ to the exact evaluation of the IRD of $V$. To yield this extension, we will choose $b$ in Algorithm 8.1 so as to yield the following bounds:$$|{\rm \tilde{det}}~V^{(q)} -{\rm det} \ V^{(q)} |< 1/2 \ , \ \ \q=2,3,\ldots ,n \ , \eqno (8.8)$$$$\parallel \tilde{V}_{\beta 1} -V_{\beta 1}\parallel <\parallelV\parallel^{-n} /2 \ , \ {\rm for \ all \ binary \ strings} \ \beta ,\ \ |\beta |< h \ , \eqno (8.9)$$where ${\rm \tilde{det}}~V^{(q)}$ and $\tilde{V}_{\beta 1}$ denote the approximations todet $V^{(q)}$ and $V_{\beta 1}$, respectively, computed by Algorithm 8.1, for the fixed value of $b$.Under (8.8) and (8.9), we recover the exact IRT of $V$ as follows:{\bf Algorithm 8.2.} {\bf Input:} a set $\{ {\rm \tilde{det}}\  V^{(q)}, \  q=1, \ldots, n\}$ of approximations to det $V^{(q)}$ for all $q$ and an approximate RD of an $n\times n$ matrix$V$ filled with integers, such that (8.8) and (8.9) hold.{\bf Output:} the IRD of $V$.\begin{enumerate}\item Round-off the values $\rm \tilde{det} \  V^{(q)}$ to the closest integers;output the resulting integer values det $V^{(q)}$, $q=2,3,\ldots ,n$.\item Compute the matrices $\tilde{W}_{\beta 1} =\tilde{V}_{\beta 1}$det $V^{(Q(\beta 1))}$ and round-off their entries to the closestintegers; output the resulting integer matrices $W_{\beta 1}=V_{\beta 1}$ det $V^{(Q(\beta 1))}$, for all binary strings $\beta$ oflength less than $h$.\itemOutput all the integer submatrices of $W_{\beta 1 \gamma}$ corresponding to thebinary strings that end with substrings $\gamma$ filled with the bits 0.\end{enumerate}{\it Correctness} of Algorithm 8.2 follows since, clearly, the values det $V^{(q)}$ are integers for all $q$ and since the matrices $W_{\beta 1} =V_{\beta 1}$det $V^{(Q(\beta 1))}$ are filled with integers for all $\beta$ (dueto Corollary 4.2). The {\it computational cost} of performing the algorithm is bounded by $O_A(1, n^2)$.{\bf Remark 8.1.} Instead of choosing the multipliers $m_{\beta 1}$based on Corollary 4.2, one may follow the more straightforwardrecipe of [R95] and recursively define $m_{\beta}$ by using induction on $|\beta|$ and by writing    $m_{\beta 1}=m_{\beta}$ det $V_{\beta 0}$. Then, however, the order of $\log |m_{\beta}|$ grows from$|\beta |$ (compare our Remark 4.1) to $|\beta |^2$, and the bit-precision and thebit-complexity  of computations grow by the extra factor $n$. [The statement of Proposition 5.1 of [R95] is false, its proof relies on an erroneous claim that if a matrix $mA$ is filled with integers, then so is the matrix $m \ {\rm adj}\ A$; this claim is false, say, for $m=3$ and the matrix $A=$ diag $(1/3,\ldots ,1/3).]$\vspace{.23in}%\noindent{\bf 9. Error and Complexity of Approximating RDs and of Computing  IRDs for a Class of C.-D.D. Matrices}We recall that, according to Definition 2.1, we   write $I$ to denote the identitymatrices of appropriate sizes. We will next specify (in terms of $n$ and $\|V\|$) a choice of the input parameter $b$ of Algorithm~8.1 that will enable us to satisfy the relations (8.8) and (8.9),  where $V$ is an$n\times n$ matrix of the following class:$$V=F-mI \ , \eqno (9.1)$$$F$ is an $n\times n$ matrix filled with non-negative integers that are less than a fixed prime $p \geq n$, and $m$ is an integer satisfying$$m=10p^2n^2 \ . \eqno (9.2)$$Due to Proposition 2.3, the entries of the matrices $V$ of this class and of all matrices $V_{\alpha}$of their RDs satisfy the following rough estimate, which will suffice for our purpose:{\bf Proposition 9.1}. {\it The entries of the matrices$V_{\alpha} +mI$ lie in the real line interval with the endpoints$2-2p$ and $2p-2$, for all binary strings $\alpha$ of length atmost $h=\log n$.}Hereafter, we will assume that $n>1$. From Proposition 9.1, we obtain two corollaries.{\bf Corollary 9.2.}$\parallel V_{\alpha 0}^{-1}\parallel < 1.1/m$, {\it for allbinary strings $\alpha$ of length at most $h-1$}.{\bf Proof.} Write $F_{\alpha 0}=V_{\alpha 0} +mI$. Due toProposition 9.1, we have $\parallel F_{\alpha 0}\parallel\leq2(p-1)n \leq (p-1)n^2$. On the other hand,$$V_{\alpha 0}^{-1}={1\over m}(I-F_{\alpha 0}/m)^{-1}=\left( {1 \over m}\right)\sum_{i=0}^{\infty} (\frac{F_{\alpha 0}}{m})^i \ .$$Therefore,$$\parallel V_{\alpha 0}^{-1}\parallel\leq \left( {1\over m}\right)\sum_{i=0}^{\infty} (\frac{\parallel F_{\alpha 0}\parallel }{m})^i< \left( {1\over m}\right) \sum_{i=0}^{\infty} {1\over (10p)^i}  = {10p \over (10p -1)m} < {1.1 \over m} \ . \hfill \Box$$Hereafter, we write$$w=m+(p-1)n \ . \eqno (9.3)$${\bf Corollary 9.3.} {\itLet $|\alpha |=h$, so that $V_{\alpha}=(v_{\alpha})$ is a $1\times 1$ matrix, and let $w$ be defined by (9.3). Then we have }$$|v_{\alpha}|\leq w \ .$${\bf Proposition 9.4.} {\it The requirement (8.8)is satisfied if the values $v_{\alpha}$, for all binary strings $\alpha$of length $h$, are approximated within an error bound }$$\delta \leq  w^{1-n}/(2n+2) \ . \eqno (9.4)$${\bf Proof.} By the virtue of Proposition 5.2, det $V^{(q)}$is the product of exactly $q$ values $v_{\alpha}$, for $\alpha$denoting binary strings of length $h$. Under the assumptions ofProposition 9.4, the errors of computing det $V^{(q)}$ are maximizedwhere $q=n$, $v_{\alpha} =w$, and the approximations to $v_{\alpha}$equal $w+\delta$ for all $\alpha$. Then, det $V=w^n$ isapproximated by $(w+\delta )^n$, with an approximationerror$$E=(w+\delta )^n -w^n= w^n ((1+(\delta /w))^n -1)=w^{n-1}\delta \sum_{i=1}^n (\delta /w)^{i-1} \left(\begin{array}{c} n \\ i \end{array}\right).$$We have $\displaystyle { \left(\begin{array}{c} n \\ 1 \end{array}\right)= n}$,  $\displaystyle{ \left(\begin{array}{c} n \\ i \end{array}\right)< 2^n}$ for all $i$, and $\delta /w<1$. Therefore, \noindent$E < (n+(\delta /w)(n-1)2^n )\delta w^{n-1}$.(9.2)--(9.4) imply that $$w>(n-1)2^n \delta \ ,$$and we may rewrite our bound on $E$ as follows:$$E<(n+1)\delta w^{n-1} \ .$$Substitute (9.4) and obtain that $E< 1/2$.\hfill $\Box${\bf Proposition 9.5.} {\itUnder (9.1)--(9.3), both requirements (8.8) and (8.9)are satisfied if the entries of the matrices $V_{\beta 1}$ for allbinary strings $\beta$ of length less than $h$ are approximatedby Algorithm 8.1 within an error bound}$$\sigma =0.5/w^n  \ . \eqno (9.5)$${\bf Proof.} From (9.1), we deduce that$$\parallel V\parallel\leq m+(p-1)(n-1)<w=m+(p-1)n  \ .$$Therefore, we will satisfy (8.9) ifwe approximate the entries of the matrices $V_{\beta 1}$ within theerror bound (9.5).On the other hand, $$\sigma =0.5/w^n <w^{1-n}/(2n+2)$$[compare (9.3)], and the values $v_{\alpha}$ for all binarystrings $\alpha$ of length $h$ are among the entries of the matrices$V_{\beta 1} $ for $|\beta | =h-1$. Therefore, the assumptions of Proposition~9.4 and, consequently, the requirement (8.8) are satisfied too. \hfill $\Box$ In Appendix A we will extend Proposition 9.5 as follows:{\bf Proposition 9.6.} {\it Under a choice of $b=O(n\log (np))$,  the bounds (8.8) and (8.9) can besatisfied  in all applicationsof Newton's iteration (7.5) within Algorithm 8.1, which is in turnapplied to approximation of the RD of a matrix $V$ satisfying (9.1) and (9.2).}So far, we have assumed infinite precision of computing the RDand IRD by means of Algorithms 8.1 and 8.2. In Appendix B, we will show thatthis is not necessary for obtaining the result of Proposition 9.6, that is, we will prove{\bf Proposition 9.7.} {\itThe estimates of Proposition 9.6 hold even if the computations byAlgorithms 8.1 and 8.2 are performed with a precision of $\tilde{b}$  bits, for some $\tilde{b}=O(n\log (np))$.}By combining the results of Corollary 5.3,Propositions 9.6, and 9.7 and by using the B-principle, we obtain{\bf Corollary 9.8}. {\itAlgorithm 8.2 computes the IRD of a matrix $V$ satisfying (9.1) and(9.2) at the cost $O_A((\log n)^2 \log (n\log p),n^{\omega}/\log n)$for $\omega$ of (6.2), (6.3); furthermore, $\tilde{b}$-bit precision suffices inthese computations for $\tilde{b}$ of order $n\log (np)$.} The computational cost estimates of Corollary 9.8 can be alternatively deduced from Proposition 6.1, since Proposition 9.1 implies that  $\| V\| \leq w$, $v \leq w$ for $v$ of (6.4), $V$ of (9.1), (9.2),  and $w$ of (9.3). On the other hand,  we only state the estimates of Corollary~9.8  for the sake of completeness of our analysis; they will be a little improved in the next section.\vspace{.23in}%\noindent\begin{center}{\bf 10. Pipelined Computation of the IRD}\end{center}Our next goal is a modification of Algorithm 8.1, which will enable us toimprove the complexity bounds of Corollary 9.8. The idea is to applythe techniques of pipelining along the lines of [PR91] (where suchtechniques have been specialized as stream contraction for computing theRD of a matrix over some semirings). Roughly speaking, we will involve thematrices $X_{\alpha ,j}$ for smaller $j$, $j < g(\alpha )$(computed well before the matrices $X_{\alpha ,g(\alpha )}$ becomeavailable) into the computation of $\tilde{V}_{\alpha 1}$ [compare (8.2), (8.6), (8.7)].More precisely, we will fix some natural $g$ and an $n \times n$ matrix $V_{\Lambda ,g}$, write $$V_{\Lambda ,i}=V_{\Lambda ,g} \    {\rm for} \   i=0, 1, \ldots, g-1 \ , \eqno (10.1)$$ where $\Lambda$ denotes the empty string, and change theequations (8.2)--(8.7) as follows: $$V_{\alpha ,i} =\left(\begin{array}{cc}B_{\alpha ,i} & C_{\alpha ,i} \\E_{\alpha ,i} & G_{\alpha ,i} \end{array}\right) \  {\rm for} \ i=0,1, \ldots, |\alpha|+g \ , \eqno (10.2)$$$$V_{\alpha ,i+1}=V_{\alpha ,i}\ , {\rm for}\ i=0,1, \ldots, |\alpha|-1 \ , \eqno (10.3) $$ $$V_{\alpha 0,i+1}=B_{\alpha ,i} \ , \ \ \ \V_{\alpha 1,i+1} =G_{\alpha ,i} -E_{\alpha ,i}Y_{\alpha ,i+1} C_{\alpha ,i} \ , {\rm for} \  i=|\alpha|, \ldots, |\alpha|+g \ , \eqno (10.4)$$where $\alpha$ ranges among all binary strings of length at most $h$  in (10.2)--(10.3) and at most $h-1$ in (10.4) and where for $i=|\alpha |,|\alpha |+1,\ldots ,|\alpha |+g$, we have $$Y_{\alpha ,i+1}=-X_{\alpha ,i,4} \ , \eqno (10.5)$$$$X_{\alpha ,i,0}=Y_{\alpha ,i} \ , \ \ {\rm for } \ i>|\alpha | \ ,\ \ \ X_{\alpha ,|\alpha |,0}=-D^{-1}(B_{\alpha ,|\alpha |} ) \ ,\eqno (10.6)$$$$X_{\alpha ,i,j}=X_{\alpha ,i,j-1} (2I+B_{\alpha ,i}X_{\alpha ,i,j-1} ) \ , \ \ \ \ \j=1,2,3,4 \eqno (10.7)$$[compare (8.3), (8.6), and (8.7)].  %% $\Lambda$ denoting the empty string, %%and $\alpha$ ranging among all binary %%strings of length less than $h$.Then, we proceed as follows:{\bf Algorithm 10.1}, {\it stream contraction for approximating the RD}.{\bf Input:} natural $g$, $h$, and $n=2^h$, and an $n \times n$ matrix $V_{\Lambda, g}$.{\bf Output:} the matrices $V_{\alpha ,|\alpha |+g}$ for allnonempty binary strings $\alpha$ of length at most $h$, satisfyingthe equations (10.1)--(10.7).{\bf Computations:}{\bf Stage 0.} Apply (10.1) to define the matrices $V_{\Lambda, i}\ , i=0, \ldots, g-1$. %%Compute $g$ of (10.2); then, starting %%with the matrix $V_{\Lambda ,0} =V$, %%inductively in $|\alpha |$ compute %%the matrices $V_{\alpha ,|\alpha |}$ %%satisfying (10.1)--(10.7), for all binary %%strings $\alpha$ of length at most $h$.{\bf Stage s}, ${\bf s=1,\ldots ,g+h}$. Starting with the matrices $V_{\alpha ,|\alpha |+s-1}$ available at the end of Stage $s-1$ for all $\alpha$ of length less than$\min (h, s)$, define the matrices $V_{\alpha 0,|\alpha |+s}$ by (10.4) and compute [concurrently in all binary strings $\alpha$ of length less than $\min(h,s)$] the matrices $V_{\alpha 1,|\alpha |+s}$ by using (10.4)--(10.7).Output the matrices $V_{\alpha ,|\alpha |+g}$ for all $\alpha$.Hereafter, for a given $n \times n$ matrix $V$, an integer $m$ of (9.2), and a positive $b$ [compare (7.2), (8.5), and (9.1)], we will assume that $g$  and $V_{\Lambda, g}$ in Algorithm~10.1 are defined as follows: $$c=\log \tilde{m} \ , \  \  \tilde{m}=\frac{m}{(n-1)(p-1)n}\ , \   \   g=1+\lceil \log ((b+2\log n)/c) \rceil \ , \eqno (10.8)$$$$V_{\Lambda ,g} =V+40 \mu m I \ , \   \   \    \mu \leq 2^{-b}/(40mn^3 ) \ . \eqno (10.9)$$Under these assumptions, in Appendix C, we will deduce{\bf Proposition 10.1.} {\it Let $p \geq n$ in (9.2). Then, the following bounds hold for all matrices$V_{\alpha ,|\alpha |+g} $ output by Algorithm 10.1:}$$\parallel V_{\alpha ,|\alpha |+g}-V_{\alpha} \parallel \leq 2^{-b} \ .$$Let us estimate the {\it computational cost} of performing Algorithm 10.1.For any $s$, Stage s is reduced essentially to 10 steps ofmultiplication of at most $n/k$ pairs of $k\times k$ matrices, for $k=2^{\ell}$, $\ell =1,2,\ldots ,h$. All these pairs of $k\times k$ matricesfor $k=2^{\ell}$ and for all $\ell$ are pairwise multipliedconcurrently, at the overall cost bounded by $O_A(\log n,n^{\omega})$, $2\leq \omega <2.376$ [compare (6.2) and observe that$\displaystyle{ \sum_{\ell =1}^h 2^{\ell} (n/2^{\ell} )^{\omega}=O(n^{\omega } )}$ for $\omega > 1$]. Summarizing these bounds for allStages s, $s=0,1,\ldots ,g+h$, we obtain{\bf Proposition 10.2.} {\itAlgorithm 10.1 supports approximating the RD of a c.-d.d. $n\times n$matrix $V$ of (9.1)  within the error norm bound $2^{-b}$, at theoverall cost $O_A((\log n)(\log n+g), n^{\omega })$for $g$ of (10.8) and $\omega $ of (6.2).}By combining Algorithms 10.1 and 8.2, summarizing the estimates for the computational costof their performance, and extending the round-off error analysis ofAppendix B, we obtain{\bf Corollary 10.3}. {\itThe IRD of a c.-d.d. $n\times n$ matrix $V$ of (9.1) can be exactlycomputed at the computational cost$O_A((\log n) \log (np), n^{\omega})$, for$\omega$of (6.2), $2\leq\omega <2.376$; moreover, this computation can beperformed by only involving operations with $\tilde{b}$-bit precisionnumbers for $\tilde{b}=O(n\log (np))$.}\vspace{.23in}\begin{center}%\noindent{\bf 11. Computing Modulo a Fixed Prime of the ERD of an IntegerMatrix}\end{center}Our next goal is probabilistic extension of Corollary 10.3 from the class of matrices $V$ of (9.1) to the class of all strongly nonsingular integer matrices $A$. In this section, we will compute the IRD and even the ERD of $A$ modulo a fixed prime $p$; in the next section we will shift to the IRD of $A$. Let $A=(a_{ij})$ be a strongly nonsingular $n\times n$ matrix filledwith integers $a_{ij}$. Then, by virtue of Proposition 3.2, there exists the RD of $A$. Let $p$ be a fixed prime, let $0 \leq f_{ij}=a_{ij}$ mod $p <p$, for all $i,j$, and let$$F=(f_{ij})=A  \ {\rm mod } \  p \ . \eqno (11.1)$$(Here and hereafter, we assume that $0\leq a$ mod $p<p$ for anyinteger $a$.)Our next goal is the computation modulo $p$ of the ERD of the matrix $F$as an auxiliary stage of computing the ERD of $A$. At first, we shouldexamine if there exists the ERD modulo $p$ of $F$.{\bf Proposition 11.1.} {\itLet $\rho >2$ be a fixed scalar, let $A$ be a strongly nonsingular $n\times n$ integer matrix, %filled with integers, where $n>1$, $\|A\|>1$, and let $p$ be a prime chosenrandomly (under the uniform probability distribution) in theinterval $J=\{p: \ n^{\rho-1}\log\parallel A\parallel <p<n^\rho \log\parallel A\parallel\}$. Then, $p \geq n$, and  the matrix F of (11.1) is stronglynonsingular modulo p with a probability at least $1-P_{\rho, n}$, for $P_{\rho, n} < Cn^{2-\rho}(1+(\rho \log n))$ and  for some positive constant C.}{\bf Proof.} It is known (see Fact 3.2.4 on p. 239 of [BP94]) that thereare at least$(n^\rho \log \parallel A \parallel )/(C \rho \log n+ \log \log  \parallel A\parallel )$ distinct primes in the interval $J$ (for some positive constant$C$). By combining this result with Proposition~2.7, one may fix $k$ and estimate the ratio between the numbers of the primes $p$ in $J$ for which (det $A^{(k)}) \bmod p=0$ and of all primes in $J$. The probability of having (det $A^{(k)}) \bmod p =0$ is given by this ratio, which is estimated in Corollary 3.2.1 of [BP94], page 240. By applying this corollary for $b=1$,$N=(\log\parallel A^{(k)}\parallel) /\log\log\parallel A^{(k)} \parallel$,we obtain that (det $A^{(k)}) \bmod p=0$ with a probability at most$P_k=Ck^{1-\rho}(1+(\rho\log k))$, for any fixed $k$, $ \ 1\leq k\leq n$.Therefore, (det $A^{(k)})\bmod p\neq 0$ for all $k$, $ \ k=1,\ldots ,n$, with a probability at least $1-P_{\rho, n}$, $$P_{\rho, n} = \sum_{k=1}^n P_k = C\sum_{k=1}^n k^{1-\rho}(1+\rho\log k)<Cn^{2-\rho}(1+\rho\log n) \ .$$\hfill $\Box$Now, we will assume that a prime $p$ has been chosen in  the interval $J$ of Proposition 11.1 and thatthe matrix $F$ of (11.1) is strongly nonsingular modulo $p$ and,therefore, possesses its ERD modulo $p$.The next algorithm computes modulo $p$ such an ERD,  representing each auxiliary or output rational value as a pair of its numerator and denominator given as two integers reduced modulo $p$. (This enable us to avoid the costly stage of computing integer reciprocals modulo $p$.) {\bf Algorithm 11.1,} {\itcomputing the ERD modulo a fixed prime.}{\bf Input:} a prime $p$ and a pair of strongly nonsingular$n\times n$ matrices $A$ and $F=A$ mod $p$ filled with integers.{\bf Output:} the (common) ERD modulo $p$ of $A$ and $F$.{\bf Computations.}{\bf Stage 0.} Compute $m=10(np)^2 $ and the c.-d.d. matrix$V=F-mI$ [cf. (9.1), (9.2)].{\bf Stage 1.} Compute the IRD of $V$ by applying Algorithms 10.1and 8.2. Then, compute modulo $p$ the RD of $V$, by dividing modulo$p$ all the computed matrices $m_{\alpha}V_{\alpha}$ of the IRD by thecomputed multipliers $m_{\alpha}$ for all $\alpha$; represent the result of each divisions by a pair of an entry of $m_{\alpha} V_{\alpha}$ reduced modulo $p$ and $m_{\alpha} \bmod p$. Output the computed RD modulo $p$ of $V$, which is also the RD modulo $p$ of $F=V$ mod $p$.{\bf Stage 2.} Recall Proposition 5.2 and compute det $V$.{\bf Stage 3.} Recall from Proposition 2.8 that $|$ det$V|\leq \parallel V\parallel^n$ and apply Newton's iteration (7.5)for $B=V$ in order to compute an approximation $X$ to $V^{-1}$satisfying (7.3) for $B=V$ and for $b$ satisfying$$ 2^{-b}/m <\parallel V\parallel^{-n} /2.2 \ . \eqno (11.2)$$Then, compute the entries of the matrix $X$ det $V$ and round them offto the closest integers, which gives us adj $V$.{\bf Stage 4.} Compute the matrix $W=({\rm adj} \ V)$ mod $p-mI$.Apply Algorithms 10.1 and 8.2 to compute the IRD of $W$. Then, compute modulo $p$ the matrices $W_{\alpha 0}$ of the RD of $W$ for all binary strings$\alpha$ of length less than $h$. Output the set of matrices,   to be denoted $\{ (W_{\alpha 0} / {\rm det} \ V) \ {\rm mod} \ p\}$ and whose entries are the pairs of integers, each reduced modulo $p$ and  such that one of them is an entry of $W_{\alpha 0} \bmod p$ and another is $(\det V) \bmod p$. (This set of matrices defines  the extending set $\{ B_{\alpha 0}^{-1}$ mod $p\}$ of the RD modulo $p$ of the input matrix $F$.)To verify correctness of Algorithm 11.1, first extend Corollary 9.2 toobtain that $\parallel V^{-1}\parallel\leq 1.1/m$. Together with (11.2),this impiles that$$\parallel X \ {\rm det} \ V-{\rm adj} \ V\parallel <1/2$$for the matrix $X$ computed at Stage 3 of Algorithm 11.1. Therefore,the rounding-off at this stage correctly defines adj $V$.Furthermore, the matrices $W_{\alpha}$ mod $p$ (see Stage 4) representthe RD modulo $p$ of adj $V$. Therefore, the set $\{ (W_{\alpha } /{\rm det} \ V) $ mod $p\}$ represents the RD modulo $p$ of $V^{-1}$.To complete the correctness proof, it remains to observe thatthe set of matrices $\{ (W_{\beta 0}/ {\rm det} \  V)\ {\rm mod} \  p,|\beta |<h\}$ is nothingelse but the extending set $\{ B_{\beta 0}^{-1}$ mod $p,|\beta |<h\}$of the (common) RD modulo $p$ of the 2 matrices $V$ and $F=V$ mod $p$.This follows from the next simple result:{\bf Proposition 11.2.} {\itLet $\{ V_{\alpha }\}$ and $\{ W_{\alpha } \}$ denote the RDs of a pair of$n\times n$ matrices V and $W=V^{-1}$, respectively. Then, $V_{\alpha}^{-1} =W_{\alpha}$ for all binary strings $\alpha$ of length at most $h$.}{\bf Proof.} Compare (2.3) and (2.4) to obtain that $V_0^{-1}=W_0$,$V_1^{-1}=S^{-1}=W_1$. Recursively extend this observation to all binarystrings $\alpha$, so as to complete the proofs of both of Proposition 11.2 and,consequently, of the correctness of Algorithm 11.1.\hfill $\Box$Similarly to deducing Corollary 10.3, we estimate the complexity ofperforming Algorithm 11.1. We arrive at{\bf Proposition 11.3.} {\itThe ERD modulo a fixed prime p of an $n\times n$ matrix A filledwith integers and strongly nonsingular modulo p (that is, such that (det $A^{(k)} $) mod $p\neq 0$, for all $k$), as well as $\det A^{(k)} \bmod p$, for all $k$, can be computed at the cost $O_A((\log n)\log (np),n^{\omega})$ for $\omega$of (6.2), $2\leq \omega <2.376$; moreover, this computation can beperformed by computing with the $\tilde{b}$-bit precision operands for$\tilde{b}=O(n\log (np))$. }{\bf Remark 11.1.} One can be tempted to simplify Algorithm 11.1 andto compute modulo $p$ the extending set $\{ V_{\alpha 0}^{-1} \}$ ofthe RD of the matrix $V$ via a more straightforward application of thetechniques of sections 3--10. In particular, one may proceed by following the recipe of [R95]: first approximate the matrices $V_{\alpha 0}^{-1}$ closely enough, then multiply the approximations by appropriateinteger multipliers $M_{\alpha}$ so as to arrive at approximations(within an error norm bound by less than 1/2) to integer matrices$M_{\alpha} V_{\alpha 0}^{-1}$, and then recover the matrices$M_{\alpha} V_{\alpha 0}^{-1}$ via rounding-off and $V_{\alpha 0}^{-1}$via divisions by $M_{\alpha}$. The problem with this approach is in bounding the size of the multipliers $M_{\alpha}$. We need to have $\log|M_{\alpha}|=\tilde{O}(n)$ in order to support the bit-precision bounds ofProposition 11.2, but if we follow the cited recipe of [R95],  we would only reach the bounds of order$\tilde{O}(n^2 )$ on $\log |M_{\alpha}|$, which would imply involving extra factor $n$ in the bit-precision and the bit-complexity bounds. Here, the notation $\tilde{O}(s)$ should be read as $O(s\log^cs)$ for a constant $c$ independent of $s$.\vspace{.23in}%\noindent\begin{center}{\bf 12. $p$-adic Lifting of the ERDs and the Recovery of the Inverses, Determinants, and Ranks of Integer Matrices}\end{center}In the previous section, we have computed the ERD modulo $p$ of an integer matrix $A$. We will now compute its $p$-adic (Newton-Hensel's) lifting, that is, the ERD modulo $p^{2^g}$ of $A$ for any fixed natural $g\geq h=\log n$. We will achieve this by incorporating the known techniques [MC79] for $p$-adic lifting of matrix inverses into our Algorithm~10.1, where we will remove the equations~(10.1) and (10.2) (since we will now fix a natural $g \geq h$ as an input value) and will modify the equations (10.4)--(10.7) so as to have: \begin{enumerate}\item $X_{\alpha, i, 4}$ replaced by $X_{\alpha, i, 1}$ in (10.5); \item $-D^{-1}(B_{\alpha, |\alpha|})$ in (10.6) replaced by $(B^{-1}_{\alpha, |\alpha|}) \bmod p$; \item the range for $j$ narrowed from $j=1, 2, 3, 4$ to $j=1$ in (10.7), and \item all arithmetic operations performed modulo $p^{2^{i-|\alpha|}}$ in (10.4) and modulo $p^{2^{i+1-|\alpha|}}$ in (10.5) and (10.7), for each pair of $\alpha$ and $i$ (as in the previous section, the input and output of arithmetic operations modulo $q$, for $q = p^j$ and for all $j$, are represented  as the ratios of pairs of integers both of which have been reduced modulo $q$). %%%, that is, the computations corresponding %%to Stage $s$ of Algorithm~10.1 performed %%modulo $p^{2^s}$, $s=1, 2, \ldots, g$.\end{enumerate}{\it Correctness} of the resulting algorithm follows since the modified equation of (10.7), that is, $$X_{\alpha, i, 1}=X_{\alpha, i, 0}(2I+ B_{\alpha, i}X_{\alpha, i, 0})\bmod p^{2^{i+1}} $$ implies that $$(I+ B_{\alpha, i}X_{\alpha, i, 1}-(I+ B_{\alpha, i}X_{\alpha, i, 0})^2) \bmod p^{2^{i+1}}=0 \ , $$and therefore, %\begin{equation}\label{equ12.1}$$Y_{\alpha, i+1} = B^{-1}_{\alpha, |\alpha|} \bmod  p^{2^{i+1}}  \eqno (12.1)%%\end{equation}$$(compare [MC79] or [BP94], Fact~3.3.1 on page 244).The {\it arithmetic complexity} estimates $O_A( (g +\log n)\log n, n^\omega)$ of Proposition~10.2 are extended to the case of the resulting algorithm (hereafter referred to as {\bf Algorithm 12.1}), where $g$ denotes a fixed natural input value, $g \geq h = \log n$. We will keep assuming that $p$ is a prime fixed in the interval $J$ of Proposition~11.1, $\|A\| > 1, \ n>1$, and the matrix $F$ of (11.1) is strongly nonsingular. Furthermore, hereafter we will assume that  $$g= 1+ \lceil\log \frac{1+n\log \|A\|}{\log p}\rceil  \eqno (12.2)$$ [cf. (10.8) and Algorithm 10.1].Then, we have  $$ 2p^4 \|A\|^n > p^{2^g} \geq 2p^2 \|A\|^n > 2\|A\|^n . \eqno (12.3)$$Therefore, by the virtue of Proposition~2.7, the value $0.5 p^{2^g}$ exceeds $|\det A|$ as well as the maximum absolute value of any entry of adj~$A$. We observe that  $$q=\left\{ \begin{array}{ll} q \bmod p  & \mbox{if $q \bmod p < 0.5 q$,} \\ (q \bmod p)-p & \mbox{otherwise,}											\end{array}   \right.$$ provided that $q$ is an integer and $2|q|<p$. These observations and (12.3) enable us to recover det~$A$ from $({\rm det}~A )\bmod p^{2^g}$ and adj~$A$ from $({\rm adj}~A )\bmod p^{2^g}$. Then, we may immediately compute $A^{-1}=({\rm adj} A)/ {\rm det}~A$, since $A$ is a nonsingular matrix.Summarizing the algorithms and the complexity estimates of this and the previous sections, we arrive at {\bf Proposition 12.1. } {\it Let $A$ be a strongly nonsingular $n \times n$ matrix filled with integers. Let $n>1$, let $\|A\|>1$, and let $p$ be a prime from the interval $J$ of Proposition~11.1 for a fixed $\rho >2$. Furthermore, let the matrix $F=A \bmod p$ be strongly nonsingular too. Then, one may compute $A^{-1}$ and $\det~A^{(k)}$, $k= 1, 2, \ldots, n$, in 2 stages that amount essentially to application of Algorithms 11.1 and 12.1, respectively, and are performed at the arithmetic cost bounded by $O_A((\log n) \log (n p), n^{\omega})$, at the first stage (compare Proposition~11.3)  and  $O_A((\log n)^2,  n^{\omega})$, at the second stage, for $\omega$ of (6.2), $2\leq \omega<2.376$.} {\bf Corollary 12.2.} {\it Under the assumptions of Proposition 12.1, one may compute $A^{-1}$ and $\det~A^{(k)}$, for $k=1, 2, \ldots, n$, at arithmetic cost $O_A((\log n) \log (n \log \|A\|), n^{\omega})$, for $\omega$ of (6.2). }Let us extend Proposition~12.1 and Corollary~12.2 so as to estimate at first the bit-precision and then the Boolean complexity of the same computations. We immediately recall the bound $O(n \log(np))$ on the bit-precision required in Algorithm~11.1, that is, at the first stage of the computations of Proposition~12.1. For the second stage (that is, essentially for Algorithm~12.1), we revisit the derivation of Proposition~10.2 and recall or estimate again that this stage is essentially reduced to $g+h$ substages, for $g$ of (12.2) and for $h=\log n$, such that the cost of performing each substage is dominated by the cost of 10 steps of multiplication of at most $n/k$ pairs of $k \times k$ matrices for $k=2^\ell$ and $\ell =0, 1, \ldots, h-1$. At each of such 10 steps, all the matrix multiplications are performed concurrently, as in the case of derivation of Proposition~10.2. Furthermore, at Substage $s$ of the second stage, for $s=1,\ldots, g+h$, the matrices of the $i$-th largest size encountered at this substage are pairwise multiplied modulo $p^{2^{t_s+1-i}}$, $i=0, 1, \ldots, h_s$, where $h_s \leq \min (h,s)$, $$t_s=s\ \  {\rm if}\ \ s\leq g, \ \ t_s=g\ \  {\rm if}\ \  s\geq g \ .\eqno (12.4)  $$   The above bounds on the moduli imply some bit-precision bounds since computation modulo $\ell$ can be performed with $\lceil \log \ell\rceil$-bit precision. Furthermore, we recall the known estimates $O_B((\log k)\log \log k, k)$ for the Boolean complexity of performing an arithmetic operation modulo $2^k-1$ (see [AHU74], [BP94], [CK91], [RT90]), which can be extended to our computations whenever we perform them with $k$-bit precision. Combining the latter estimates with ones for the arithmetic cost and for the bit-precision of our computations, we bound the Boolean cost of performing Algorithm~11.1, that is,  the first stage of the computations of Proposition~12.1 (cf. Proposition 11.3), by  $O_B((\log n) \log (np) \log (n\log p)\log \log (n \log p), n^{\omega+1}\log (np))$ and the  \\ Boolean cost of performing the $s$-th stage of Algorithm~12.1  by \\$O_B((\log n)(\log (2^{t_s}\log p)\log \log (2^{t_s}\log p), n^{\omega+1}2^{t_s}\log p)$ [cf. (12.4)], $s=1, \ldots, g+h$.  (Recall that the $s$-th stage of Algorithm 12.1 is the $s$-th substage of the second stage of the computations of Proposition~12.1.)  By summarizing all these estimates, for $p$ lying in the interval $J$ of Proposition~11.1 and for $g$ satisfying  (12.2), (12.3), we estimate the Boolean complexity of our computations. To simplify the expressions for the resulting estimates, we write $$A=(a_{i,j}), \ \  a = \log \max_{i,j}|a_{i,j}| \ , \eqno (12.5) $$so that $g=O(\log (na))$, $g+h = O(\log (na))$,   $\log (2^g \log p)=O(\log (na))$ for $g$ of (12.2), $\log p =O(\log (na))$, $\log(np) = O(\log (na))$, $\log (n \log p)=O(\log (n \log a))$.Then, we rewrite our Boolean cost bounds as follows: \\$O_B((\log n)(\log (na))(\log (n\log a))\log \log (n \log a), n^{\omega+1}\log (na))$ for performing  Algorithm 11.1, $O_B((\log n)(\log (na))\log \log (n a), n^{\omega+1}a)$ for performing the $s$-th stage of Algorithm 12.1 for $s= g+1, \ldots, g+h$, and \\   $O_B((\log n)(s +\log \log (na))\log (s + \log \log (n a)), n^{\omega+1}2^s \log (na))$ for performing the $s$-th stage of Algorithm 12.1 for $s=1, \ldots, g$. By applying the B-principle, we bound the overall cost of performing the first $g=O(\log(na))$ stages of Algorithm 12.1 by \\$O_B((\log n)(\log (na))^2 \log \log (na), n^{\omega +1} a /\log (na))$, and we bound the overall cost of performing its last $h=\log n$ stages by \\$O_B((\log n)^2(\log(na))\log \log (na), n^{\omega +1}a)$. Then again, we apply the B-principle, to yield the same parallel Boolean time bound, that is, \\$O((\log n)(\log (na))^2 \log \log (na))$, in all the 3 estimates (for Algorithm 11.1, for the first $g$ stages of Algorithm 12.1, and for its last $h$ stages), which gives us the Boolean processor bounds  \\  $O((n^{\omega +1}\log(n \log a) \log \log (n \log a))/ \log \log (n a))= O(n^{\omega +1} \log (n \log  a))$,\\  $O(n^{\omega +1}a /\log (na))$, and  $O((n^{\omega +1} a \log n)/\log (na))$, for these 3 groups of computations, respectively.By using the Boolean cost bounds of Proposition 11.3 for computing $\det A^{(k)} \bmod p$, for all $k$, and by combining the cited Boolean time bound and the three latter processor bounds, we obtain {\bf Proposition 12.3.}  {\it Under the assumptions of Proposition~12.1, one may compute the inverse matrix $A^{-1} \bmod p$ and $\det A^{(k)} \bmod p$, $k=1, \ldots, n$, at the Boolean cost bounded by \\ $O_B((\log n)(\log (na))(\log (n \log a)) \log \log (n \log a), n^{\omega+1}\log (na))$  and the matrix $A^{-1}$ and $\det A^{(k)}$, $k=1, \ldots, n$, at the Boolean cost bounded by \\$O_B((\log n)(\log (na))^2\log \log (n a), n^{\omega+1} (\frac{a }{\log (n a)}+1)\log n)$  for $\omega$ of (6.2) and $a$ of (12.5).} {\bf Remark 12.1.} The estimates of Proposition 12.3 and Corollary 12.2 have been stated in terms of the upper bound $na$ on the magnitudes of the logarithms of $\det A$ and the entries of $A^{-1}$ but can be restated in terms of any other available upper bound on these magnitudes (which can be substantially smaller than $na$  for some classes of the input matrices $A$). A similar observation applies to the estimates of Theorems 1.1 and 1.2. Let us extend the complexity estimates of Corollary~12.2 and Proposition~12.3 to ones of Theorem~1.1 [note that $n \log \|A\| = O(n a)$, where  $2^a=\max_{i,j}|a_{ij}|$, according to (12.5)]. We will first relax the assumptions about strong nonsingularity of $A$. Towards this goal, we fix $\rho >2$, fix a random prime $p$ in the interval $J$ of Proposition~11.1, and extend Algorithm~11.1 so as to compute $({\rm det}~(AA^T)^{(k)}) \bmod p$, for $k=1, \ldots, n$, and $r(p)=\max\{k, \  \det~(AA^T)^{(k)} \bmod p \neq 0\}$.Let us write $r= \max\{k, \  \det~(AA^T)^{(k)} \neq 0\}$, so that $r={\rm rank} \  A$ (due to Proposition~2.10). Then, clearly, $r\geq r(p)$. Furthermore, $r=r(p)$, with a probability $1- P_{\rho, n}$, estimated in Proposition 11.1.Thus, we output $r(p)$ as rank~$A$ and arrive at the estimate of Theorem 1.1 for the randomized cost of computing rank~$A$. If $r(p)<n$, then we output $\det A =0$. Otherwise, $\det A \neq 0$, and $AA^T$ is a strongly nonsingular matrix. Then, by the virtue of Proposition~11.1, $AA^T$ is strongly nonsingular modulo $p$, with a probability $1-P_{\rho,n}$. Therefore, we may apply the results of this section so as to   compute at first $(AA^T)^{-1}$ and then $A^{-1}=A^T(AA^T)^{-1}$ and $\vec{x}= A^{-1}\vec{f}$ satisfying $A\vec{x}=\vec{f}$. [Strong nonsingularity (modulo $p$) of $AA^T$ is tested as a by-product of computing $(AA^T)^{-1}$.] We may also immediately compute $\det (AA^T) = (\det A)^2$, though this does not give us the sign of $\det A$.    Now, assume that $r(p)<n$. Then, the $r(p) \times r(p)$ leading principal submatrix $B=(AA^T)^{(r(p))}$ of $AA^T$ is nonsingular [GL89]. Let us write $AA^T=\left( \begin{array}{lr} B & C\\ D & E  \end{array} \right)$, $G=\left(  \begin{array}{lr} I & -B^{-1}C \\ O & I \end{array} \right)$,and observe that $AA^TG=\left(  \begin{array}{lr} B & O \\ D & Q \end{array} \right)$, where $Q=O$ if and only if $r(p)={\rm rank} \  A$. (Compare [KP91] and [BP94], pp. 110 and 333, though we now use symmetrization rather than randomization.) This gives us an algorithm for verification whether $r(p)=r$ (at the cost within the asymptotic cost bounds of Theorem 1.1). If $r(p)=r$, indeed, then the $n-r$ columns of the matrix $A^TG\left(  \begin{array}{c} O\\ I  \end{array} \right)$, where $I$ denotes the $(n-r) \times (n-r)$ identity matrix, give us a basis for the null-space, ${\bf N}(A)$,  of $A$ (compare Definition 2.7). We recall from fact 2.12 that if there exists a solution $\vec{x}$ to a linear system $A\vec{x}=\vec{f}$, then it  can be represented as $\vec{x}=\vec{x_0}+\vec{z}$, $\vec{x_0}$ being a fixed specific solution and $\vec{z}$ being a vector from ${\bf N}(A)$. Let $\vec{g}$ be the $r$-dimensional prefix-subvector of $\vec{f}$, made by the first $r$ components of $\vec{f}$.  Let $\vec{y}=B^{-1}\vec{g}$ be the solution to the nonsingular system $B\vec{y}=\vec{g}$. Then, a specific solution $\vec{x_0}$ to the system $A\vec{x}=\vec{f}$ is given by $\vec{x_0}=A^TG\left(  \begin{array}{c} \vec{y}\\ \vec{0} \end{array} \right)$ if the latter linear system is consistent, and we have $A\vec{x} \neq \vec{f}$ otherwise. To complete the proof of Theorem 1.1, it remains to furnish an algorithm for computing $\det A$. We will do this in the next section.\vspace{.23in}%\noindent\begin{center}{\bf 13. Some Definitions and Auxiliary Results on Computations with Structured Matrices}\end{center}We will recall some definitions and some simple and/or well-known facts (cf. [BP94], [CKL-A87], [KKM79], [P92]), partially repeating section 1.1.{\bf Definition 13.1.} $T=(t_{i,j})$ is a $k\times k$ {\it Toeplitz  matrix} if $t_{i+1, j+1} = t_{i,j}$ for $i,j = 0, 1, \ldots, k-2$ [cf. (1.1)], that is, if the entries of $T$ are invariant in their shifts in the diagonal direction. (Such a matrix is defined by its 2 columns (or rows), the first one and the last one.) A lower triangular Toeplitz matrix $L$ is defined by its first column ${\bf u}$ and is denoted $L({\bf u})$. $Z=Z_k=(z_{i,j})$ is a $k\times k$ lower triangular Toeplitz matrix with the first column $(0,1, 0, \ldots, 0)^T$, so that $z_{i+1, i}=1$ for all $i$, $i=0, 1, \ldots, k-2$, $z_{i,j}=0$ if $j\neq i+1$. {\bf Proposition 13.1.} {\it The product of a $k\times k$ Toeplitz matrix and a vector of dimension $k$ can be computed at the cost $O(\log k, k)$ [via reduction to 3 FFTs, each on $O(k)$ points, or to convolution of 2 vectors of dimension $O(k)$].}{\bf Definition 13.2.} For a $k \times k$ matrix $A$, write $F_+(A) =A-ZAZ^T$, $F_-(A)=A-Z^TAZ$. If $F(A)=GH^T$ for a pair of $k \times \ell$ matrices $G$ and $H$ and for $F=F_+$ or $F=F_-$, then the pair of $G, \ H$ is called an $F$-{\it generator} of $A$ of length $\ell$.  (Note that, in this case, the pair $H$, $G$ is an $F$-generator of $A^T$ of the same length.) The minimum length $\ell$ of an $F$-generator of $A$, for fixed $A$ and $F$, is called the $F$-{\it rank} of $A$, is denoted $r_F(A)$, and is equal to ${\rm rank}~F(A)$. A $k \times k$ matrix $A$ is called a {\it Toeplitz-like matrix} if it is given with its $F$-generator (for $F=F_+$ or $F=F_-$) having a length bounded by a constant independent of $k$. {\bf Proposition 13.2.} {\it $r_F(T) \leq 2$ if $T$ is a Toeplitz matrix, and $r_F(T) \leq 1$ if $T$ is a triangular Toeplitz matrix, for $F=F_+$ and $F=F_-$. In particular, $r_F(I)=1$.} The correlation to definition 1.2 is given by the following result. {\bf Proposition 13.3.} {\it $G, \ H$ is an $F_+$-generator (respectively, $F_-$-generator) of $A$ having a length $\ell$, $G=({\bf g_1}, \ldots, {\bf g_{\ell}})$,  $H=({\bf h_1}, \ldots, {\bf h_{\ell}})$,  if and only if \\ $A=\sum_{s=1}^{\ell} \   L({\bf g_s})~L^T({\bf h_s})$ [respectively, if and only if $A=\sum_{s=1}^{\ell} \L^T({\bf g_s})~L({\bf h_s})$].}An $F_+$-generator (respectively, $F_-$-generator) of a length $\ell$ for any matrix $A$ can be immediately transformed [at the cost $O_A(\log n, n)$ of performing $O(1)$ FFTs] into an $F_-$-generator (respectively, $F_+$-generator) of length at most $\ell +2$ for $A$.Based on the latter results, we will operate with the $F$-generators of Toeplitz-like matrices, rather than with the matrices themselves. Such a representation is storage efficient and also enables  us to use less time and/or fewer processors in Toeplitz-like computations, due to the following corollary  (cf. Propositions~13.1  and 13.3):{\bf Corollary~13.4.} {\it The product of a $k \times k$ Toeplitz-like matrix by a vector of dimension $k$ can be computed at the cost $O_A(\log k, k)$.}We also have {\bf Proposition 13.5.} {\it Given $F$-generators, $G_A, \ H_A$ of length $\ell_A$ and $G_B, \ H_B$ of length $\ell_B$ of $k\times k$ matrices $A$ and $B$, respectively, for $F=F_+$ or $F=F_-$, one may compute an $F$-generator $G_{AB}, \ H_{AB}$ of $AB$ of length at most $\ell_A + \ell_B +1$ at the cost $O_A(\log k, (\ell_A+\ell_B)^2k)$, and an $F$-generator of $A+B$ of length at most $\ell_A + \ell_B $ is immediately available cost-free.} {\bf Proposition 13.6.} {\it $r_+(A)\leq r_-(A) +2$, $r_-(A)\leq r_+(A) +2$  for any matrix $A$. If $A$ is nonsingular, then $r_+(A^{-1}) = r_-(A)$.}{\bf Proposition 13.7}(cf. Proposition~A.6 of [P92]). {\it Let an $F$-generator of a $k \times k$ matrix $A$ of length $\ell$ (for $F=F_+$ or $F=F_-$) and an upper bound $r^*$ on the $F$-rank $r_F(A)$ be given. Then, an $F$-generator of $A$ of length at most $r^*$ can be computed at the cost $O_A(r^*,~k~r^*)$.} The next result is immediately verified (compare Definition~2.4). {\bf Proposition 13.8.} {\it Let $GH^T=F_+(W)$ for a $k\times k$ matrix $W$. Then $(GH^T)^{(i)}=F_+(W^{(i)})$ for $i=1, 2, \ldots, k$; furthermore, $r_{F_+}(C) \leq r_{F_+}(W)+1$,  $r_{F_+}(E)  \leq r_{F_+}(W) +1 $, under (2.1), and  $r_{F_+}(T) \leq r_{F_+}(W) + 2$ for any submatrix $T$ of $W$.}It follows that $r_{F_+}(B) \leq r_{F_+}(W)$, under (2.1). We observe similar relations for trailing principal submatrices and the operator $F_-$. (2.4) implies that $S^{-1}$ is a trailing principal submatrix of $W^{-1}$. Therefore, $r_{F_-}(S^{-1}) \leq r_{F_-}(W^{-1})$. By applying Proposition 13.6 for $A=S$ and $A=W$, we obtain that $r_{F_+}(S) \leq r_{F_+}(W)$. {\bf Proposition 13.9.} {\it Let (2.1) and (2.2) hold, where $B$, $S$, and $W$ are nonsingular matices. Then $\max \{ r_{F_+}(B) , r_{ F_+ }(S) \} \leq r_{F_+}(W) $.}By recursively applying the latter proposition, we obtain {\bf Corollary 13.10.} {\it Let $V_{\alpha}$ be a matrix of the RD of a matrix $A$. Then, $r_{F_+}(V_{\alpha}) \leq r_{F_+}(A)$.} We will use the following result from [KS91] (cf. also Lemmas 1.5.1 and 2.13.1 of [BP94]).{\bf Proposition 13.11.} {\it Let $S$ be fixed finite set of cardinality $|S|$. Let $A$, $L$, and $U$ be $n \times n$ matrices, rank $A=r$, $L$ and $U^T$ being unit lower triangular Toeplitz matrices, defined by the $2n -2$ entries of their first columns. Let these entries be chosen from $S$ at random, independently of each other, under the uniform probability distribution on $S$. Then the matrix $(UAL)^{(r)}$ is strongly nonsingular with a probability at least $1-(r+1)r/|S|$.}  By using the randomization of Proposition 13.11, we may compute $\det A=\det (UAL)$, since $UAL$ is strongly nonsingular, with a probability at least $1-(n+1)n/|S|$ if $A$ is nonsingular. By letting $|S|=n^4$, say, and by applying Propositions 11.1 and 12.3 to the matrix $UAL$, we arrive at the desired algorithm for $\det A$, which enables us to complete the proof of Theorem 1.1. \hfill $\Box$Let us next extend the algorithm supporting Proposition~13.7.We will achieve this by applying the {\it SVD truncation of displacement generator}.{\bf Algorithm 13.1} [P93]. {\bf Input:} $F=F_+$ or $F=F_-$, an $F$-generator $G$, $H$ of a $k \times k$ matrix $A$ of length $r$, and a natural $r'<r$.{\bf Output:} an $F$-generator $G'$, $H'$ of a $k \times k$ matrix $A$ of length at most $r'$ such that $$ \|A' - A\|_2 \leq (1+2(r_{F_+}(A)-r')k) \min_Y \|Y-A\|_2  \ , \eqno (13.1)$$where the minimum is over all $k\times k$ matrices $Y$ of $F$-rank at most $r'$.{\bf Computations.}{\bf Stage 1.} Compute the singular value decomposition (SDV) of the matrix $GH^T=F(A)$, that is, compute a pair $U$ and $V$ of unitary $k \times k$ matrices and a  $r \times r$ diagonal matrix $\Sigma = {\rm diag}(\sigma_1, \ldots, \sigma_r)$ for positive $\sigma_1, \ldots, \sigma_r$ satisfying $$ GH^T = F(A) = U \Sigma  V^T. $$ {\bf Stage 2.} Compute and output an $F$-generator $G'$, $H'$ of $A'$ of length at most $r'$  as follows: $$G'= U \Sigma _{r'}, \ \  H' = VI_{r', r} , $$where $\Sigma _{r'} = {\rm diag}(\sigma_1, \ldots, \sigma_{r'}, 0, \ldots, 0)$ and  $I_{r',r}= {\rm diag}(1, \ldots, 1, 0, \ldots, 0)$ are $r\times r$ matrices of rank $r'$. {\it Correctness} of this algorithm has been proved in [P93]. The {\it computational cost} of its performance [where the SVD is computed approximately and the requirement (13.1) is weakened by factor 2 (say)] is bounded by $O_A(\log n, n/\log n)$ for $r=O(1)$, [P93a].{\bf Remark 13.1.} If $r' \geq r_F(A)$, then (13.1) implies that $\|A'-A \|_2 = \min_Y \|Y-A\|_2=0$, and then Algorithm~13.1 is an equivalent alternative to the algorithm supporting Proposition~13.7, except that the latter algorithm is rational (it can be performed with no errors over the rationals), whereas Algorithm~13.1 has a nonrational, though numerically stable stage of computing the SVD.This suggests that the algorithm supporting Proposition ~13.7 should be applied in Algorithm~12.1,  at the $p$-adic lifting stage, whereas Algorithm~13.1 is abetter candidate to use in numerical applications of Algorithm~11.1, when it is performed with rounding-off.\vspace{.23in}%\noindent \begin{center}{\bf 14. Improvement of the Algorithms for the ERD, IRD, Inverse, Determinant, and Rank in the Toeplitz-like Case}\end{center}Let us re-examine the computation of the ERD and IRD of a strongly nonsingular $n \times n$ matrix $A$ filled with integers in the case where $A$ is a Toeplitz-like matrix given with its $F$-generator $G$, $H$ of length $r=r_{F_+}(A)$ (compare Proposition 13.7).We recall that $r_{F_+}(V_{\alpha}) \leq r$ for all matrices $V_{\alpha}$ of the RD of $A$ (compare Corollary~13.10), and we will apply either the algorithm supporting Proposition~13.7 or Algorithm~13.1 in order to decrease (to a level at most $r$) the length of the computed $F$-generators of these matrices, in all cases where this length exceeds $r$. Likewise, we will apply Proposition~13.8 and obtain that the computation of $V_{\alpha 1, i+1}$ according to (10.4) only involves matrices whose $F$-ranks are bounded from above by $3r+r_{F_+}(Y_{\alpha, i+1})+6$. According to our analysis, the matrix $Y_{\alpha, i+1}$ approximates $B^{-1}_{\alpha, i}$, and since $r_{F_+}(B_{\alpha, i}) \leq r$, we have $r_{F_-}(B_{\alpha, i}^{-1}) \leq r, \ \ r_{F_+}(B_{\alpha, i}^{-1}) \leq r+2$ (compare Proposition~13.6). We will apply Proposition~13.7 and/or  Algorithm~13.1 in order to compute an $F_+$-generator of length at most $r+2$ for a matrix $Y'_{\alpha, i+1}$ approximating $Y_{\alpha, i+1}$ and, therefore, also $V_{\alpha, i}^{-1}$. [The approximation of  $V_{\alpha, i}^{-1}$ by $Y'_{\alpha, i+1}$ slightly deteriorates versus the approximation by $Y_{\alpha, i+1}$, but since the matrix $V_{\alpha, i}$ is c.-d.d. and since $d(V_{\alpha, i})$ is small (compare Proposition~9.1),  we more than compensate for such a deterioration by performing an extra Newton step in (10.7).] Then, all matrices involved in the computation of the ERD and the IRD of $A$ will be represented by their $F_+$-generators of length  $O(r)$. Similar argument is applied to the computation of the $p$-adic lifting of the ERD of $A$, except that this argument is simplified since  (12.1) and Proposition~13.6 together imply that $$r_{F_+}(Y_{\alpha, i+1} \bmod p^{2^{i+1}}) \leq r_{F_-}(Y_{\alpha, i+1} \bmod p^{2^{i+1}}) +2 = $$ $$r_{F_+}(B^{-1}_{\alpha, |\alpha|} \bmod p^{2^{i+1}}) + 2 \leq r+2, $$ and we keep the length of the associated $F_+$-generators bounded by means of applying the rational algorithm that supports Proposition~13.7, instead of applying  Algorithm~13.1. In fact, we may also apply yet an alternative technique for bounding the length of an $F$-generator of $Y_{\alpha, i+1}$; such a technique relies on using distinct operators $F$, such as $F^+(A)=AZ-ZA$ (see [BP94], pp. 189).Finally, it is easily verified (cf. [P,b]) that the computation (of section 12) of a basis for the null-space of $A$  also involves only matrices represented by their $F_+$-generators of length $O(r)$, for a matrix $A$ given with its $F_+$-generator of length $r$.Let us now turn to estimating the computational cost, in the case of Toeplitz or Toeplitz-like input. The algorithms supporting Proposition~13.7 and Algorithm~13.1 accompany matrix multiplications that we count in the proof of Proposition~10.2. Thus, these algorithms are applied $O(2^i \log (n \log \|A\|))$ times to $(n/2^i) \times (n/2^i)$ matrices, $i=1, \ldots, h$, and therefore, are performed at the overall cost \\$O_A((\log n) \log (n\log \|A\|), n\log n)$ if $r=O(1)$. The same estimates apply to the overall cost of performing all Toeplitz-like matrix multiplications in Algorithms~11.1 and 12.1 [compare the equations (10.4) and (10.7)] because we operate with $F_+$-generators of Toeplitz-like matrices involved in these multiplications and apply Propositions~13.3, 13.5, and Corollary~13.4. This enables us to decrease the arithmetic processor bound of Corollary~12.3 from $n^{\omega}$ to $n \log n$, since concurrent multiplications of $O(2^i)$ pairs of $(n/2^i) \times (n/2^i)$ Toeplitz-like matrices, for $i=1, \ldots, h$, $h=\log n$, are performed at the overall cost bounded by $O_A(\log n, n \log n)$ and since we apply the latter bounds at most $g+h=O(\log (n \log \|A\|))$ times. The estimated overall cost of all applications of Toeplitz-like matrix multiplications,  of Algorithm~13.1, and of the algorithm  supporting Proposition~13.7 clearly dominates the cost of performing all other operations.Summarizing, we obtain the following result: {\bf Proposition 14.1.} {\it If the $n \times n$ input  Toeplitz-like matrix $A$ is filled with integers, then one may modify the computation of its ERD and IRD according to the algorithms of sections 8--12, so as to perform these computations at the overall cost $O_A((\log n) \log (n \log \|A\|), n\log n)$.} To obtain a similar extension of the Boolean complexity bounds  of Proposition~12.3 and Theorem~1.1, let us examine the precision of computations by our algorithm simplified  in the Toeplitz-like case. We recall that the fast and processor efficient parallel computation of convolution of a pair of $k$-dimensional vectors can be performed at the cost $O_A(\log k, k \log \log k)$ over any field or ring of constants ([CK91] or [BP94],  Theorem~1.7.1 on p. 56). Since effective Toeplitz-like computations can be ultimately reduced to vector convolution (Propositions~13.1 and 13.5), we may perform them with the precision of $O(2^i \log p)$ bits when we seek the output modulo $p^{2^i}$ (for any $i$ in the range from 1 to $g$). Therefore, we may perform the $p$-adic lifting stage of our algorithm with the same bit-precision as in the case of an arbitrary integer input matrix. When we approximate the ERD of an input Toeplitz-like matrix, we will effectively reduce the computations to performing FFTs (see Propositions~13.1 and 13.3) and will recall Corollary~3.4.1 of pp. 255--256 of [BP94], which shows a numerically stable implementation of FFT. We also recall that the known algorithms for the computation of the SVD of a matrix are numerically stable (see [GL89], [P93]). From these observations, we deduce that we may perform the computations with the same bit-precision (up to a constant factor independent of $n$) in both cases of our original Algorithm~11.1 for an arbitrary $n \times n$ input matrix and of its Toeplitz-like modification. Since in the latter case, we use by factor $n^{\omega -1}/\log n$ fewer arithmetic processors, we will also use by factor $n^{\omega -1}/\log n$ fewer Boolean processors in this case, thus replacing $n^{\omega}$ of (6.2) by $n \log n$ in Corollary~12.3.This enables us to extend Theorem~1.1 so as to arrive at Theorem~1.2. {\bf Remark 14.1.} In [R95], an algorithm of [BA80] is used in order to decrease the length of an $F$-generator of a matrix $A$ to the level $r = {\rm rank}\ F(A)$. Unlike our SVD truncation Algorithm~13.1 and our algorithm supporting Proposition~13.7, the algorithm of [BA80] only works if $(F(A))^{(r)}$, the $r \times r$ l.p.s. of $F(A)$, is nonsingular. Moreover, the error estimates of Proposition~7.4 of [R95] require that the matrix $(F(A))^{(r)}$ be well-conditioned. The proof of the latter proposition in [R95] relies on the assumption that $(F(A))^{(r)}$ is well-conditioned (and therefore, is nonsingular) as long as $A$ is a c.-d.d. matrix. Furthermore, on p. 753 (after Lemma 2.1) of [R95], it is proposed to decrease the length of an $F$-generator of $A$ to the level $r={\rm rank} \ F(A)$ based on the claim that $(F(A))^{(r)}$ is a nonsingular matrix for any matrix $A$. The cited assumption is, however, wrong [consider, for instance, the matrix $A={\rm diag}(I,-I)$ where $I$ is the $k \times k$ identity matrix, $k>1$, $r=2$, $(F_+(A))^{(r)} = \left( \begin{array}{lr} 1 & 0 \\ 0 & 0 \end{array} \right)$], and the cited claim is wrong as well [consider, for instance,  $A= \left( \begin{array}{lr} 0 & 0 \\ 0 & 1 \end{array} \right)$, $F=F_+$]. Furthermore, nontrivial techniques are required in order to fix these errors, and no such techniques are available in [R95], so that both algorithm and claimed complexity bounds of [R95] fall apart. {\bf Remark 14.2.} In [R95], the author has overlooked the fact of the growth of $F$-ranks of the approximations to  the inverse $A^{-1}$ of a Toeplitz-like matrix $A$ computed in the process of Newton's iteration. (Was he misled by the fact that the $F$-rank of $A^{-1}$ is $O(1)$ and has wrongly extrapolated this bound to the $F$-ranks of the approximations?) As a result, the algorithm proposed in [R95] contains no means (unlike, say, our Algorithm~13.1) in  order to prevent the length of the $F$-generators from growing to the level $n$ in less than $\log n$ Newton's steps, even if the weak statement of Lemma 2.4 (b) of [R95] is strengthened by the known record estimates.  (In fact, had we relied  on the excessively large upper estimates of the part (b) of Lemma 2.4 of [R95], we would have had the length of the $F$-generators increasing  to the level $n$ already in $O(\log  \log n)$ Newton's steps.) This growth immediately implies the growth by extra factor $n^{\omega -1}/\log n$, for $\omega$ of (6.2), $\omega \geq 2$, of (both arithmetic and Boolean) processor complexity bounds, versus ones claimed in [R95].{\bf 15. Discussion.} Our paper leaves an open question of theoretical importance, whether the level of our parallel complexity estimates of Theorem 1.2 for Toeplitz and Toeplitz-like computations can be reached by means of purely algebraic approach, with using no rounding-off to the closest integers. This question is also of practical interest because the algorithms presented in this paper require  to use the precision of order $\log |\det A|$, which generally means the order $n \log A$, even if we only need the output with a much lower precision. Historically, similar open problem had arisen for computations with general integer matrices, after appearance of [P85], [P87], and in that case (for general integer matrices), the subsequent works of [KP91], [KP92], [P91], [P92] have given us an alternative algebraic solution with no use of rounding-off. Will this be eventually done also in the Toeplitz-like case or at least in the Toeplitz case?\begin{thebibliography}{92} \footnotesize \bibitem[AGr88]{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.	 \bibitem[AHU74]{AHU74}	A. V. Aho, J. E. Hopcroft, J. D. Ullman, {\it The Design and Analysis of Computer Algorithms,} Addison-Wesley, Reading, Massachusetts, 1974.\bibitem[B68]{B68} E. H. Bareiss, Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination, {\it Math. Comput.,} 22, 565--578, 1968.	 \bibitem[BA80]{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.\bibitem[Be68]{Be68}		E. R. Berlekamp, {\it Algebraic Coding Theory,} McGraw-Hill, New York, 1968.\bibitem[Be84]{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.%%\bibitem[BG74]{bg74} A. Ben-Israel, %%T. N. E. Grevil, {\it Generalized Inverses: %%Theory and Applications}, Wiley, New York, 1974.\bibitem [BGH82]{bgh82}  A. Borodin, J. von zur Gathen, J. Hopcroft,  Fast Parallel Matrix and GCD Computation,  {\it  Information and Control}, 52, 3, 241--256,  1982.  \bibitem[BGY80]{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.\bibitem[BK87]{BK87}		A. Bruckstein, T. Kailath, An Inverse Scattering Framework for Several Problems in Signal Processing, {\it IEEE ASSP Magazine,} 6-20, January 1987.\bibitem[BL80]{BL80}  D. Bini, G. Lotti, Stability of Fast Algorithms for Matrix Multiplication,  {\it Numer. Math.,} 36, 63--72,  1980. \bibitem[BP91]{BP91}  D. Bini, V. Y. Pan,  Parallel Complexity of Tridiagonal Symmetric Eigenvalue Problem, {\it Proc. 2nd Ann. ACM-SIAM Symposium on DiscreteAlgorithms,} 384--393, ACM Press, New York, and SIAM Publications, Philadelphia, Pennsylvania, 1991.\bibitem[BP94]{bp94}D. Bini, V.Y. Pan, {\it Polynomial and Matrix Computations}, vol. 1, Birkhauser, Boston, 1994.\bibitem[BSS89]{BSS89}  L. Blum, M. Shub, S. Smale,  On a Theory of Computations Over theReal Numbers: NP-Completeness,  Recursive Functions and UniversalMachines,  {\it  Bull. of Amer. Math. Soc.,} 21,  1--46, 1989. 	  \bibitem[BT71]{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. \bibitem[BT90]{BT90}		M. Ben-Or, P. Tiwari, Simple Algorithm for Approximating All Roots of a Polynomial with Real Roots, {\it J. of Complexity,} 6, 417--442, 1990.\bibitem[Ch85]{ch84}  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.  \bibitem[CK91]{CK91} D. G. Cantor, E. Kaltofen, On Fast Multiplication of Polynomials over Arbitrary Rings, {\it Acta Inf.,} 28, 7, 697--701, 1991.  \bibitem[CKL-A87]{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.\bibitem[Cs76]{cs76}  L. Csanky,  Fast Parallel Matrix Inversion Algorithm,  {\it  SIAM J. Comput.,}   5, 4,  618--623,  1976.\bibitem[dH87]{dH87}  F.R. deHoog,  On the Solution of Toeplitz Systems,  {\it  Linear Algebra  Appl.,} 88/89,  123--138,  1987. \bibitem[E67]{E67} J. Edmonds, Systems of Distinct Representatives and Linear Algebra, {\it J. Res. Nat. Bur. Stand.,} 71B(4), 241--245, 1967. \bibitem[EG88]{EG88}  D. Eppstein, Z. Galil,  Parallel Algorithmic Techniques for Combinatorial Computation,  {\it  Annual Review of Computer Science,} 3,  233--283,  1988.	  \bibitem [G84]{g84}  J. von zur Gathen,  Parallel Algorithms for Algebraic Problems,{\it  SIAM J.  Comput.},  13, 4,  802--824,  1984.   \bibitem [G86]{g86}  J. von zur Gathen,  Parallel Arithmetic Computations: A Survey, {\it Proc. Math. Foundation of Computer Science, Lecture Notes in Computer Science},  233, 93--112, Springer, Berlin, 1986.\bibitem[GKO94]{GKO94}		I. Gohberg, T. Kailath, V. Olshevsky, Fast Gaussian Elimination with Partial Pivoting for Matrices with Displacement Structure, preprint, {\it Information Science Lab., Stanford University,} Stanford, California, 1994.\bibitem[GL89]{GL89}  G.H. Golub, C.F. Van Loan, {\it Matrix Computations,} Johns Hopkins Univ. Press,  Baltimore,  Maryland,  1989.	  \bibitem[H94]{He94}		G. Heinig, Inversion of Generalized Cauchy Matrices and Other Classes of Structured Matrices, {\it IMA Volume: Linear Algebra for Signal Processing in Math. and Its Applications,} 69, 95--114, Springer, 1994.\bibitem[J92]{J92}			J. J\`{a} J\`{a}, {\it An Introduction to Parallel Algorithms,} Addison-Wesley, Reading, Massachussetts, 1992.\bibitem[K74]{K74}			T. Kailath, A View of Three Decades of Linear Filtering Theory, {IEEE Trans. Inform. Theory,} 20, 146--181, 1974.\bibitem[K87]{K87}  T. Kailath,  Signal Processing Applications of Some Moment Problems, {\it  Proc. AMS Symp. in Applied Math.,} 37,  71--100,  1987.  \bibitem[KKM79]{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.\bibitem[KP91]{KP91}		E. Kaltofen, V. Y. Pan, Processor Efficient Parallel Solution of Linear Systems over an Abstract Field, {\it Proc. 3rd Ann. ACM Symp. on Parallel Algorithms and Architectures,} 180-191, ACM Press, New York, 1991.\bibitem[KP92]{KP92}		E. Kaltofen, V. Y. Pan, Processor-Efficient Parallel Solution of Linear Systems II.  The Positive Characteristic and Singular Cases, {\it Proc. of 33rd Ann. IEEE Symp. on Foundations of Computer Science,} 714--723, IEEE Computer Society Press, 1992.\bibitem [KP94]{KP94}E. Kaltofen, V.Y. Pan,Parallel Solution of Toeplitz and Toeplitz-like Linear Systems over Fields of 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.\bibitem [KR90]{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.\bibitem[KS91] {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.  \bibitem[KVM78]{KVM78}		T. Kailath, A. Viera, M. Morf, Inverses of Toeplitz Operators, Innovations, and Orthogonal Polynomials, {\it SIAM Rev.,} 20, 1, 106--119, 1978.\bibitem[L-AK84]{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.\bibitem[L-AKC84]{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.\bibitem[Le92]{Le92}		F. T. Leighton, {\it Introduction to Parallel Algorithms and Architectures: Arrays, Trees \& Hypercubes,} Morgan Kaufmann Publishers, San Mateo, California, 1992.	 \bibitem[Ma75]{Ma75}		J. Makhoul, Linear Prediciton: A Tutorial Review, {\it Proc. IEEE,} 63, 4, 561--580, 1975.\bibitem[Mc79]{Mc79}  R.T. Moenck, J.H. Carter,  Approximate Algorithms to Derive Exact Solutions to Systems of Linear Equations,  {\it  Proc. EUROSAM,  Lecture Notes in Computer Science,}  72,  63--73,  Springer, Berlin, 1979. \bibitem[Morf80]{Morf80} M. Morf,  Doubling Algorithms for Toeplitz and Related Equations, {\it  Proc. IEEE Internat. Conf. on ASSP}, 954--959, IEEE Computer Society Press, 1980.\bibitem[Mu81]{Mu81} B.R. Musicus, Levinson and Fast Choleski Algorithms for Toeplitzand Almost Toeplitz Matrices, \rm Internal Report, {\it Lab. ofElectronics, M.I.T.,}  Cambridge, Massachusetts,  1981.\bibitem[P85]{P85}			V. Y. Pan, Fast and Efficient Parallel Algorithms for the Exact Inversion of Integer Matrices, {\it Proc. 5th Ann. Conf. on Foundations of Software Technology and Theoretical Computer Science, Lecture Notes in Computer Science,} 206, 504--521, Springer Verlag, 1985. \bibitem[P87]{P87}		V.Y. Pan, Complexity of Parallel Matrix Computations, {\it Theoretical Computer Science,} 54, 65--85, 1987.\bibitem[P90]{P90}			V. Y. Pan, Computations with Dense Structured Matrices, {\it Math. of Computations,} 55, 191, 179--190, 1990.\bibitem[P91]{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, {\it NATO ASI Series, Series F: Computer and Systems Sciences,} 77, 27-56, Springer, Berlin, 1991.\bibitem[P92]{P92}		V. Y. Pan, Parametrization of Newton's Iteration for Computations with Structured Matrices and Applications, {\it Computers and Mathematics (with Applications),} 24, 3, 61--75, 1992.\bibitem[P92a]{P92a}			V. Y. Pan, Complexity of Computations with Matrices and Polynomials, {\it SIAM Review,} 34, 2, 225--262, 1992.\bibitem[P92b]{P92b}		V. Y. Pan, Parallel Solution of Toeplitz-like Linear Systems, {\it J. of Complexity,} 8, 1--21, 1992.\bibitem[P93]{P93}		V. Y. Pan, Decreasing the Displacement Rank of a Matrix, {\it SIAM J. on Matrix Analysis,} 14, 1, 118--121, 1993.\bibitem[P93a]{P93a}			V. Y. Pan, Concurrent Iterative Algorithm for Toeplitz-like Linear Systems, {\it IEEE Trans. on Parallel and Distributed Systems,} 4, 5, 592--600, 1993.\bibitem[P93b]{P93b}		V. Y. Pan, {\it Parallel Solution of Sparse Linear and Path Systems, in Synthesis of Parallel Algorithms} (J.H. Reif editor), Chapter 14, pp. 621--678.  Morgan Kaufmann publishers, San Mateo, California, 1993.\bibitem[P95]{P95}			V. Y. Pan, Optimal (up to Polylog Factors) Sequential and Parallel Algorithms for Approximating Complex Polynomial Zeros, {\it Proc. 27th Ann. ACM Symposium on Theory of Computing,} 741--750, ACM Press, New York, 1995.\bibitem[P96]{P96}			V. Y. Pan, A New Approach to Parallel Computation of Polynomial GCD and to Related Parallel Computations over Abstract Fields, {\it Proc. Seventh Annual ACM-SIAM Symposium on Discrete Algorithms,} 518--527, ACM Press, New York, and SIAM Publications, Philadelphia, PA, 1996.\bibitem[P96,a]{P96a}			V. Y. Pan, Optimal and Nearly Optimal Algorithms for Approximating Polynomial Zeros, {\it Computers \& Math. (with Applications),} to appear in 1996.\bibitem[P96,b]{P96b}V. Y. Pan, Parallel Computation of Polynomial GCD and Some Related Parallel Computations over Abstract Fields, {\it Theor. Computer Science}, to appear in 1996.\bibitem[P96,c]{P96c} V. Y. Pan,  Effective Parallel Computations with Toeplitz and  Toeplitz-like Matrices Filled with Integers, {\it LAM: Proc. of Summer Seminar on Math. of Numerical Analysis: Real Number Algorithms},  (Edited by J. Renegar, M. Shub, S. Smale) 1996. \bibitem[PP92]{PP92}		V. Y. Pan, F. P. Preparata, Supereffective Slow-down of Parallel Computations, {\it Proc. 4th Ann. ACM Symp. on Parallel Algorithms and Architectures,} 402-409, ACM Press, 1992.\bibitem[PP95]{PP95}		V. Y. Pan, F. P. Preparata, Work-Preserving Speed-up of Parallel Matrix Computations, {\it SIAM J. on Computing,} 24, 4, 811--821, 1995. \bibitem[PR91]{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.\bibitem[PR93]{PR93}		V. Y. Pan, J. Reif, Fast and Efficient Parallel Solution of Sparse Linear Systems, {\it SIAM J. on Computing,} 22, 6, 1227-1250, 1993.\bibitem[PSLT93]{PSLT93}		V. Y. Pan, A. Sadikou, E. Landowne, O. Tiga, A New Approach to Fast Polynomial Interpolation and Multipoint Evaluation, {\it Computers \& Math. (with Applications),} 25, 9, 25--30, 1993.\bibitem[Q94]{Q94} M.J. Quinn, {\it Parallel Computing: Theory and Practice,} McGraw-Hill, New York, 1994. \bibitem[R93]{R93}  J. 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.\bibitem[R95]{R95}			J. 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, NY, 1995.\bibitem[RT90]{RT90}  J. Reif, S.R. Tate,  Optimal Size Division Circuits,  {\it SIAM J.Comput.,} 19, 5, 912--925, 1990.\bibitem[St69]{St69}  V. Strassen,  Gaussian Elimination is Not Optimal,  {\it  Numer. Math.,}  13,  354--356,  1969.\bibitem [Ste94]{Ste94} 		W. F. Stewart, {\it Introduction to the Numerical Solution of Markov Chains,} Princeton University Press, New Jersey, 1994.\bibitem[W65]{W65} J. H. Wilkinson, {\it The Algebraic Eigenvalue Problem,} Clarendon Press,  Oxford,  1965.\end{thebibliography}\vspace{.23in}%\newpage\begin{center}{\bf APPENDIX A. Proof of Proposition 9.6}\end{center}Consider a path in the tree $T$ from the root $V$ to a leaf$V_{\alpha}=(v_{\alpha} )$, $|\alpha |=h$. Algorithm 8.1 follows such a path, by recursively going from matrices$V_{\beta}$ to $V_{\beta 0}$ and $V_{\beta 1}$. For given $V_{\beta}$and, consequently, $V_{\beta 0}$, the algorithm approximates the matrices$V_{\beta 0}^{-1}$, within error norm bounds $2^{-b}\parallelV_{\beta 0}^{-1}\parallel$, and then extends such an approximationto approximating $V_{\beta 1}$. The errors of the computed approximation to$V_{\beta}$ are extended and accumulated in computation of all thedescendands of $V_{\beta}$ along the paths in $T$. We need toestimate the resulting overall errors in all the output matrices alongall such paths, assuming that $b$ is large though bounded by being of order$O(n\log (np))$.Hereafter, we will write $\tilde{M}$ to denote the approximationsto matrices $M$ computed by Algorithm 8.1, for $M$ denoting$V_{\alpha}$ (for any binary string $\alpha$), a submatrix of$V_{\alpha}$, or any other auxiliary matrix involved. We will also write$$\Delta (M)=\tilde{M}-M \ . \eqno (A.1)$$We will first estimate $\parallel\Delta (V_{\alpha 1} )\parallel$ in termsof $\parallel\Delta (V_{\alpha })\parallel$. For convenience, we write$W=V_{\alpha}$, $S=V_{\alpha 1}$, and recall (2.1), (2.2). FromProposition 9.1 and Corollarly 9.2, we obtain that$$\max \{\parallel B\parallel ,\parallel C\parallel ,\parallel E\parallel , \parallel G\parallel \} \leq \parallel W\parallel \leq w \ , \eqno (A.2)$$$$\parallel B^{-1}\parallel\leq 1.1/m \ . \eqno (A.3)$$We also write [cf. (A.1)]:$$\tilde{W}=\left(\begin{array}{cc}\tilde{B} & \tilde{C} \\\tilde{E} & \tilde{G} \end{array}\right) , \ \ \\Delta (W) =\left( \begin{array}{cc}\Delta (B) & \Delta (C) \\\Delta (E) & \Delta (G) \end{array}\right) \ ,$$$$\tilde{S}=\tilde{G}-\tilde{E} L\tilde{C} \ ,$$where $L$ denotes the computed approximation to $\tilde{B}^{-1}$.Then, clearly,$$\max\{ \parallel \tilde{B}\parallel ,\parallel \tilde{C}\parallel ,\parallel \tilde{E}\parallel ,\parallel \tilde{G}\parallel \}\leq \parallel \tilde{W}\parallel , \eqno (A.4)$$$$\max \{ \parallel \Delta (B)\parallel ,\parallel \Delta (C)\parallel ,\parallel \Delta (E)\parallel ,\parallel \Delta (G)\parallel \}\leq \parallel \Delta (W)\parallel \ . \eqno (A.5)$$We will assume that the errors of approximations obtained viaAlgorithm 8.1 are sufficiently small so that the following inequalitieshold:$$\parallel L-\tilde{B}^{-1}\parallel\leq \nu\parallel\tilde{B}^{-1}\parallel , \ \ \ \ \nu\leq 1/(4000m^2 ) , \eqno (A.6)$$$$\parallel\tilde{B}^{-1}\parallel < 1.11/m , \ \ \\parallel L\parallel\leq (1+\nu )\parallel\tilde{B}^{-1}\parallel < 1.11 (1+\nu )/m . \eqno (A.7)$${\bf Remark A.1.} The first inequality is reconciled with (7.3) if$\nu\leq 2^{-b}$.{\bf Proposition A.1.} {\itSuppose that the inequalities (A.2)--(A.7) hold and that$$40 \nu m\leq\parallel \Delta (W)\parallel < 0.01  \ w , \eqno (A.8)$$where $W=V_{\alpha}$ and $\alpha$ is a binary string of length lessthan $h$. Then,}{\rm a)} $\ \\parallel \Delta (B)\parallel =\parallel\Delta (V_{\alpha 0} )\parallel \leq \parallel\Delta (V_{\alpha })\parallel =\parallel\Delta (W)\parallel \ , ${\rm b)} $ \ \\parallel\Delta (S)\parallel =\parallel\Delta (V_{\alpha 1})\parallel \leq 5 \parallel \Delta (V_{\alpha } )\parallel =5\parallel\Delta (W)\parallel$.{\bf Proof.} Part a) of the proposition follows immediately since$V_{\alpha 0}$ is a submatrix of $V_{\alpha}$. To deduce part b), we needthe following bound [implied by (A.2) and (A.8)]:$$\parallel\tilde{W}\parallel < 1.01 \ w \ , \eqno (A.9)$$as well as the next simple proposition and corollary:{\bf Proposition A.2.} {\itFor any 2 pairs of $k\times k$ matrices, $X$ and $\tilde{X}$ and $Y$ and $\tilde{Y}$, we have}a) $ \ \ \Delta (X\pm Y)=\Delta (X) \pm \Delta (Y)$,b) $ \ \ \Delta (XY)=\Delta (X)Y+X\Delta (Y)+\Delta (X)\Delta (Y)$,\noindent{\it and if $X$ and $\tilde{X}$ are nonsingular matrices, then also}c) $ \ \ \Delta (X^{-1})=X^{-1}\Delta (X)\tilde{X}^{-1}=-\tilde{X}^{-1}\Delta (X)X^{-1}$.{\bf Corollary A.3.} {\itUnder the assumptions of Proposition A.2, we have}a) $ \ \ \parallel\Delta (X\pm Y )\parallel \leq \parallel\Delta (X) \parallel +\parallel \Delta (Y) \parallel$,b) $ \ \ \parallel\Delta (XY)\parallel\leq\parallel\Delta (X)\parallel \\parallel Y\parallel +\parallel\Delta (Y) \parallel \\parallel X\parallel +\parallel\Delta (X) \parallel \\parallel \Delta (Y) \parallel$,c) $ \ \ \parallel\Delta (X^{-1} )\parallel \leq \parallel X^{-1}\parallel \ \parallel \tilde{X}^{-1}\parallel \\parallel \Delta (X)\parallel$.Next, we apply part c) of Corollary A.3 for $X=B$ and $\Delta (X^{-1})=\tilde{B}^{-1}-B^{-1}$ and obtain that$$\parallel\Delta (B^{-1} )\parallel\leq \parallel B^{-1}\parallel \\parallel \tilde{B}^{-1}\parallel \ \parallel\Delta (B)\parallel \ .$$Substitute (A.3), (A.4), (A.5), and (A.7) into the latter bound and obtain that $$\parallel\Delta (B^{-1})\parallel < 1.221 \parallel\Delta (W)\parallel /m^2 \ .$$Extend this bound by applying (A.6), (A.7), and (A.8) and arriveat the relations$$\parallel L-B^{-1}\parallel\leq \parallel \Delta (B^{-1})\parallel+\parallel L-\tilde{B}^{-1}\parallel\leq \parallel\Delta (W)\parallel(1.221 + (1.11)/40 )/ m^2$$$$\leq 1.3 \parallel\Delta (W)\parallel /m^2 \ . \eqno (A.10)$$Apply part b) of Corollary A.3 for $X=E$, $Y=L$ and deduce that$$\parallel\Delta (EL)\parallel\leq \parallel \Delta (E) \parallel \\parallel L \parallel + \parallel L-B^{-1}\parallel (\parallel E\parallel+\parallel\Delta (E) \parallel ) \ .$$Recall from (9.2) and (9.3) that $w/m \leq 1.02$. By combining thetwo latter inequalities with our bounds on $\parallel L\parallel$,$\parallel L-B^{-1}\parallel$, $\parallel E\parallel$, and$\parallel\Delta (E)\parallel$ [see (A.4)--(A.10)], obtain that$$\parallel\Delta (EL)\parallel\leq \parallel\Delta (W)\parallel \left({1.11\over m} (1+\nu ) +\left( {1.3\over m^2} \right) 1.01 w\right) \leq \frac{2.5}{m} \parallel\Delta (W)\parallel  \ . \eqno (A.11)$$Then again, we apply part (b) of Corollary A.3, this time for$X=EL$, $Y=C$, and obtain that$$\parallel\Delta (ELC)\parallel \leq \parallel\Delta (EL)\parallel \\parallel \tilde{C}\parallel +\parallel\Delta (C)\parallel \\parallel EL\parallel \ .$$Substitute our previous estimates (A.4)--(A.9), (A.11) into the latterinequality and deduce that$$\parallel\Delta (ELC)\parallel\leq\parallel\Delta (W)\parallel\left( \left({2.5\over m}\right) 1.01 w+{1.11\over m} (1+\nu )w\right) \leq$$$$3.7\parallel\Delta (W) \parallel w/m  \ ,$$and since $w/m\leq 1.02$, we have$$\parallel\Delta (ELC)\parallel\leq 4\parallel\Delta (W)\parallel \ .$$From (A.5), we have $\parallel\Delta (G)\parallel\leq\parallel \Delta(W)\parallel$. Therefore, by applying part a) of Corollary A.3 for$X =G$, $Y=ELC$, we deduce that$$\parallel\Delta (S)\parallel\leq\parallel\Delta (ELC)\parallel+\parallel \Delta (G)\parallel\leq 5\parallel\Delta (W)\parallel \ ,$$which proves Proposition A.1.\hfill $\Box$To satisfy the assumptions of Proposition A.1, we first set$W=V$, $\tilde{W}=V+40\nu \ mI$, so that$\parallel\Delta (W)\parallel =40\nu m$, thus satisfying (A.8).[We could have set $W=\tilde{W}=V$, but this would have violated (A.8).]We also write $L=-X_3$, where $X_0=D^{-1} (\tilde{W})$ and where $X_1$, $X_2$, and $X_3$ are defined by (7.5). Then, (7.6) impliesthat $$\parallel X_i +\tilde{B}^{-1}\parallel\leq\parallel\tilde{B}^{-1}\parallel /m^{2^i},\  \ i=0,1,2,3,$$ and consequently,$$\parallel L-\tilde{B}^{-1}\parallel\leq \nu\parallel\tilde{B}^{-1}\parallel, \ \ for \ \ \nu\leq 1/m^8 < 1/(4000 m^2 n^{12} ),$$ thus satisfying(A.6). The remaining assumption(A.7) of Proposition A.1 is also easily verified. By applyingProposition A.1 recursively, we recursively extend its assumptions to$W=\tilde{V}_0$, $W=\tilde{V}_1 +40 \nu m I$, and $W=\tilde{V}_{\alpha}+40 \nu m I$ for any binary string $\alpha$of length at most $h$.In particular, we obtain that$$\parallel\Delta (W)\parallel\leq 40 \nu m \sum_{i=0}^{|\alpha |} 5^i < 40 \nu m n^3\ ,  $$ for $W =\tilde{V}_{\alpha} +40 \nu m I$ and for all $\alpha$ (with $|\alpha |\leq h=\log n)$, which we easily reconcilewith (A.6)--(A.8). Recall Remark A.1, substitute $\nu =2^{-b}$, andobtain that$$\parallel\Delta (V_{\alpha})\parallel <(40 mn^3 )2^{-b} <2^{2-b}w^2n \ , \eqno (A.12) $$for all $\alpha$.Let us choose $b$ of order $n\log (np)$ so as to satisfy the bound$$b\geq 3 + \log n +(n+2)\log w,$$ which is compatible with the choice of $b=\log(1/\nu)$ and with (A.6). Substitute this bound on $b$ into the preceding upper bound on $\parallel\Delta(V_{\alpha})\parallel$ and obtain that$$\parallel\Delta (V_{\alpha})\parallel <0.5 \ w^{-n} \ ,$$ for all $\alpha$, satisfying the requirement of Proposition 9.5. Thiscompletes the proof of Proposition 9.6.\hfill $\Box$\vspace{.23in}%\noindent\begin{center}{\bf APPENDIX B. Proof of Proposition 9.7}\end{center}Let us first assume that the computations of Newton's step (7.5) areperformed with infinite precision, but the transition from $W=V_{\alpha}$to $S=V_{\alpha 1}$ for all binary strings $\alpha$, $|\alpha |<h$,is performed by means of computations with rounding-off to $\tilde{b}$-bitprecision. The latter transition involves two multiplications andone subtraction of $k\times k$ matrices for $k=2^{h-|\alpha |}$ [compare(2.2)]. By applying the techiques of backward error analysis [W65],[BL80], we bound the norm of the matrix $\epsilon (S)$ of the errors of the approximation to $S$caused by rounding-off, and we obtain that $$\parallel\epsilon (S)\parallel\leq n^{O(1)} 2^{-\tilde{b}} \parallelW\parallel (1+\parallel L\parallel \ \parallel W\parallel ) \ ,$$for $L$ denoting the computed approximation to $B^{-1}$ [cf. (A.6), (A.7)]. By applying the bounds (A.2), (A.7), (9.2),(9.3), and Proposition 9.1, we obtain that$$\parallel\epsilon (S)\parallel\leq m^{O(1)} 2^{-\tilde{b} } \ .$$By choosing $\tilde{b}$ of order $n\log (np)$, we make $\parallel\epsilon (S)\parallel$ less than $2\parallel\Delta (W) \parallel$ for$W=V_{\alpha}$ and for all $\alpha$. This is less than one half of the upper bound that we have in part b) of Proposition A.1. Combining both of these boundsgives us cumulative upper bound $7\parallel\Delta (W)\parallel$,reflecting the impact of the above round-off errors. This enables us topreserve the validity of the bound (A.12) (since $\sum_{i=0}^h 7^i < n^3$ for $h=\log n$) and, consequently, the validity of the entire proof of Proposition 9.6.It remains to estimate the impact of rounding-off to $\tilde{b}$-bitprecision in performing Newton's steps (7.5). Then again, we deal with 2 matrixmultiplications [if we ignore the simple addition step $2I+(WX_{i-1})$]. Byapplying backward error analysis again, we estimate that$$\parallel\epsilon (X_i )\parallel\leq n^{O(1)} 2^{-\tilde{b}}\parallel X_{i-1}\parallel (2+\parallel W\parallel \\parallel X_{i-1}\parallel ) \ , \eqno (B.1)$$where $\epsilon (X_i )$ denotes the matrix of the errors ofapproximation of $X_i$ due to rounding-off in performing the iteration (7.5). We deduce from (8.3) and Proposition 9.1 that $\| I+BX_0\| \leq \frac{1}{\tilde{m}}$, for $\tilde{m} = m/((n-1)(p-1)n) > 10p $ and for $X_0 =T/m$. Then, we obtain from  (7.6) that$$\parallel I+BX_{i-1}\parallel\leq 1/\tilde{m}^{2^{i-1}} \ ,$$and therefore,$$\parallel X_{i-1} +B^{-1}\parallel\leq \parallel B^{-1}\parallel /\tilde{m}^{2^{i-1}} \ .$$Consequently,$$\parallel X_{i-1}\parallel \leq \parallel B^{-1}\parallel\left( 1+{1\over \tilde{m}^{2^{i-1} } } \right) \ ,$$for $i=1,2,\ldots$. Substitute (A.3) and arrive at the inequality$\parallel X_{i-1}\parallel\leq  1.1(1+1/\tilde{m})/m < 1.21/m$, for $i\geq 1$.Substitute this bounds on $\parallel X_{i-1}\parallel$ and the bound$\parallel W \parallel \leq w$ [see (A.2)] into (B.1) and obtain that$$\parallel\epsilon (X_i )\parallel\leq (np)^{O(1)} 2^{-\tilde{b}} \ ,\ \ {\rm for \ all} \ i \ .$$By choosing a sufficiently large $\tilde{b}$, though of order $n\log(np)$, we easily ensure that, say,$$\parallel\epsilon (X_i )\parallel <2^{-b}/\parallel B\parallel \ .$$Therefore,$$\parallel I+B(X_i +\epsilon (X_i ))\parallel\leq \parallel I+BX_i \parallel +\parallel B\parallel \\parallel\epsilon (X_i )\parallel \leq$$$$\parallel I+BX_{i-1}\parallel^2 + 2^{-b} \ .$$Since Newton's iteration (7.5) stops if $\parallel I+BX_{i-1} \parallel \leq 2^{-b}$, we may assume that $\parallel I+BX_{i-1} \parallel >2^{-b}$, so that rounding-off may at worst change (7.6) intothe bound$$\parallel I+B(X_i +\epsilon (X_i ))\parallel\leq 2\parallel I+BX_i \parallel\leq 2^{\sigma} \parallel I+BX_0 \parallel^{2^i} < 2(2\parallel I+BX_0\parallel)^{2^i}, $$where $\sigma =\sum_{s=0}^i 2^s$.Since $\parallel I+BX_0 \parallel \leq 1/m$, the impact of the rounding-off is more thancompensated by a single extra step (7.5). Thus, doubling the number ofNewton's steps (7.5) more than compensates the effect of rounding-off,and this completes the proof of Proposition 9.7.\hfill $\Box$\vspace{.23in}%\newpage\begin{center}{\bf APPENDIX C: Proof of Proposition 10.1}\end{center}We will largely follow the line of Appendix A. Hereafter, we write$\Delta (V_{\alpha ,i} )=V_{\alpha ,i} -V_{\alpha}$ for all$\alpha$ and $i,i\geq |\alpha |$ [compare (A.1)].{\bf Proposition C.1.} {\itFor every pair $\alpha$ and $i,i\geq |\alpha |$, $|\alpha | < h$, wehave}a) $ \ \ \parallel \Delta (V_{\alpha 0,i+1})\parallel\leq\parallel\Delta (V_{\alpha ,i})\parallel$,b) $ \ \ \parallel \Delta (V_{\alpha 1,i+1})\parallel\leq 5\parallel\Delta (V_{\alpha ,i})\parallel$.The proof is similar to one for Proposition A.1. The validity of thenecessary assumptions [similar to (A.2)--(A.7)] is also verifiedinductively, but in this case, we need double induction: at first in$|\alpha |$, for the fixed $i=|\alpha |$, and then in $i$, for eachfixed $\alpha$. The induction in $i$ is almost trivial, because, clearly,the transition from $V_{\alpha ,i}$ to $V_{\alpha ,i+1}$ only decreasesthe approximation error norms $\parallel\Delta (V_{\alpha ,i})\parallel$,due to the stage of Newton's iteration (10.7) [compare (7.6)].Let us briefly comment on the induction in $|\alpha |$ for$i=|\alpha |$, which goes similarly to the induction used in Appendix A. From Proposition9.1 and the second equation of (10.6), we obtain that$$\parallel I +B_{\alpha ,|\alpha |} \ X_{\alpha ,|\alpha |,0} \parallel \leq (n-1)(p-1)n/m < 0.1/p \ ,$$for all $\alpha$. By extending (7.6) to the case of Newton's step(10.7), we obtain that$$\parallel I +B_{\alpha ,|\alpha |} \ X_{\alpha ,|\alpha |,4}\parallel\leq 1/(10p)^{16} \ ,  \ \ {\rm for \ all} \ \alpha \ .$$Since, by assumption, $p>n$, it is not hard to extend the latter bound,by using the inductive assumption of Proposition C.1,so as to arrive at anappropriate extension of (A.6) (compare theanalysis in Appendix A). By using (10.9) we extend (A.8) for$\alpha =\Lambda$, and then we extend (A.8) to all $\alpha$ by applyingProposition C.1 inductively in $|\alpha |$.Now, we obtain from (10.9) that$$\parallel V_{\Lambda ,g} -V\parallel=40\nu m \leq 2^{-b}/n^3 \ .$$By applying Proposition C.1 inductively in $|\alpha |$, we obtain that$$\parallel V_{\alpha ,|\alpha |+g} -V_{\alpha}\parallel \leq40\nu m \sum_{i=0}^{|\alpha|}5^i \parallel V_{\Lambda ,g} -V\parallel <$$$$\left(5^{|\alpha| +1}/(4n^3) \right)  2^{-b} < 2^{-b} \ ,$$for all binary strings $\alpha$ of length at most $h=\log n$. Thisproves Proposition 10.1.\hfill $\Box$\vspace{.23in}%\newpage\begin{center}{\bf APPENDIX D: Outline of Some Techniques of [P85], [P87].}\end{center}The basic idea of [P85], [P87] is to shift from the integer input matrix $V$ to the auxiliary matrix $F=V+ apI$, where $I$ is the identity matrix, $p$ is a prime, $a$ is a natural number, and $p$ and $a$ are sufficiently large. For large $p$ and $a$, the matrix $-F$ is strongly diagonally dominant, and its inverse, $-F^{-1}$, is well approximated by $X_0=-I/(ap)$. By applying Newton's iteration (7.5) for $B=F$, %%$$ X_{i+1}=X_i(2I +FX_i),~~i=%%0, 1, \ldots, \eqno (D.3) $$ one may very rapidly improve this approximation and very closely approximate $-F^{-1}$. Since $F$ is an integer matrix, $\det F$ and the entries of adj~$F = (\det F)F^{-1}$ are integers, which can be recovered by rounding-off their approximations within absolute errors less than 1/2. The remaining ingredient is the approximation of $\det F$. In [P85], [P87], this is achieved as a by-product of computing $\det(xI-F)$ the characteristic polynomial of $F$. In the present paper, we employ a more routine approach, via recursive block factorization of $F$ (cf. [St69], [AHU74], [BA80], [Morf80], [PR91], [P93], [PR93]). Having computed $F^{-1}=V^{-1} \bmod p$ and $\det F = \det V \bmod p$, we apply $p$-adic (Newton-Hensel's) lifting and rapidly arrive at $\det V$, $V^{-1}$ and adj~$V$ (modulo a high power of $p$) and then recover $\det V$, adj~$V$ and $V^{-1}$.   	   		 \end{document}