From:	IN%"Victor.Pan@sophia.inria.fr"  "Victor Pan"  3-APR-1997 11:23:14.41To:	IN%"vpan@korrigan.inria.fr"CC:	Subj:	determinant file (with Y.Y> and C.S.)Return-path: <Victor.Pan@sophia.inria.fr>Received: from vice1.lehman.cuny.edu by lcvax.lehman.cuny.edu (PMDF #12809) id <01IH9MGF0XDC000DE1@lcvax.lehman.cuny.edu>; Thu, 3 Apr 1997 11:22 ESTReceived: from korrigan.inria.fr by vice1.lehman.cuny.edu (5.51/SMI-3.0DEV3) id AA06313; Thu, 3 Apr 97 10:57:13 ESTReceived: by korrigan.inria.fr (8.8.5/8.6.12) id SAA25288; Thu, 3 Apr 1997 18:14:36 +0200 (MET DST)Date: Thu, 03 Apr 1997 18:14:36 +0200From: Victor Pan <Victor.Pan@sophia.inria.fr>Subject: determinant file (with Y.Y> and C.S.)To: vpan@korrigan.inria.frMessage-id: <199704031614.SAA25288@korrigan.inria.fr>Content-type: text/plain; charset=us-asciiX-Mailer: exmh version 1.6.7 5/3/96Mime-Version: 1.0------- Forwarded MessageReturn-Path: YING_AND_STEVE@msn.comReceived: from sophia.inria.fr by korrigan.inria.fr (8.8.5/8.6.12) with ESMTP id EAA13827 for <vpan@korrigan.inria.fr>; Fri, 28 Mar 1997 04:37:14 +0100 (MET)Received: from upsmot02.msn.com by sophia.inria.fr (8.8.5/8.6.12) with SMTP id EAA20958 for <Victor.Pan@sophia.inria.fr>; Fri, 28 Mar 1997 04:35:31 +0100 (MET)Received: from upmajb04.msn.com ([204.95.110.81]) by upsmot02.msn.com (8.6.8.1/Configuration 4) with SMTP id TAA07537 for <Victor.Pan@sophia.inria.fr>; Thu, 27 Mar 1997 19:31:48 -0800Date: Fri, 28 Mar 97 01:26:46 UTFrom: "Steve Yu" <YING_AND_STEVE@msn.com>Message-Id: <UPMAIL09.199703280137230626@msn.com>To: "Victor Pan" <Victor.Pan@sophia.inria.fr>Subject: RE: file\documentstyle [12pt]{article}\topmargin -.10in\oddsidemargin .05in\evensidemargin .05in\textheight 8.50in\textwidth 6.25in\newtheorem {algorithm}{\bf Algorithm}[section]\newtheorem {problem}{\bf Problem}[section]\newtheorem {theorem}{\bf Theorem}[section]\newtheorem {proposition}{\bf Proposition}[section]\newtheorem {lemma}{\bf Lemma}[section]\newtheorem {corollary}{\bf Corollary}[section]\newtheorem {fact}{\bf Fact}[section]\newtheorem {remark}{\bf Remark}[section]\newtheorem {definition}{\bf Definition}[section]\newtheorem {claim}{\bf Claim}[section]\renewcommand{\baselinestretch}{1.5}\def\theequation{ \arabic{section}.\arabic{equation} }\begin{document}\addtocounter{page}{-1}   %%%%%{-11}%%\addtocounter{page}{+6} %%%{+6}\title{\bf \Large ALGEBRAIC AND NUMERICAL TECHNIQUES FORTHE COMPUTATION OF MATRIX DETERMINANTS }\noindent\author{\small  Victor  Y.  Pan\thanks{Supported by NSF Grants  CCR 9020690 and CCR 9625344 andby PSC CUNY Awards 665301, 666327 and 667340.} \\{\small Department of Mathematics and Computer Science}\\{\small Lehman College, City University of New York,Bronx, NY 10468}\\{\small Internet:  VPAN@LCVAX.LEHMAN.CUNY.EDU}\\{\small Yanqiang Yu$^*$  }\\{\small Colin Stewart$^*$ }\\{\small Ph.D. Program in Computer Science}\\{\small Graduate Center of City University of New York, New York, NY 10038}}%%%\thispagestyle{empty}\date{}%%%\pagestyle{myheadings}%%%\markboth{V. Y. Pan}{Algebraic and Numerical Techniques for%%%the Determinant Computation }\maketitle\noindent\begin{abstract} { \small\baselineskip=9ptWe review, modify and combine together several numerical and algebraictechniques in order to compute the determinant of amatrix or the sign of such a determinant.The resulting algorithms enable us to obtain the solution by usinga lower precision of computations and a relatively few arithmeticoperations. The problem has importantapplications to computational geometry. }\vspace{.3in}\noindent{\bf Key Words:} evaluation of the determinant,the sign of the determinant, matrix singularitytest, modular (residue) arithmetic, rounding error analysis.\vspace{.3in}\noindent{\bf 1991 Math. Subject Classification:} 68U05, 65Y25, 52B55, 14Q99, 68Q40,15A15\end{abstract}\newpage\noindent\section{\bf Introduction.}\label{sec1}\setcounter{equation}{0}\subsection{\bf The subject and some background.}We study the classical problems of the computationof the determinant of a matrix or testing if thedeterminant vanishes, that is, if the matrix issingular. These problems have a long history (see,for instance, [Mu], [Ma], [Du], [Fo], [R52], [E67],[B68], [P84], [P88] and [BP]) and have recently received anew major motivation, due to their importantapplications to geometric computations, suchas computation of convex hulls and Voronoi diagramsand testing if the line intervals of a given familyhave a non-empty common intersection. In suchapplications, one needs {\it sign} or {\it singularity}tests, that is, one needs either to test%if det $A = 0$ or, more generally,if det $A > 0$, det $A = 0$, or det $A < 0$, foran $n \times n$ matrix $A$, or just to testwhether $\det A=0$ or not. In one group of theseapplications, $n$ is {\it relatively small} [C], [ABDPY],ranging from 2 to 10, but mostly staying below 5, and theentries of $A$ represent high precision approximations toreal data. In this class of applications, computations areusually performed with the 64-bit double precision (whereeither floating point arithmetic withabout 50 bits allocated to the mantissas is used orfixed point representation with the double precision of 64bits is used), though [ABDPY] also proposes a singleprecision approach. In another major group of geometricapplications [DLa], [ES], [AF], [FR], [MA], oneneeds to know orientation of a high-dimensionalpolyhedron or of a high-dimensional algebraic variety (say, for theproblems of convex optimization in statistical physicsand chemistry). Then the range for $n$ is on a muchhigher level, say, from 100 to 500, whereas the entriesof the input matrix $A$ can be represented with a lowerprecision of 5 to 10 bits. In both of these groups ofapplications, the matrix $A$ can be either fixed or updateddynamically, so that every time only one of its columns(or one of its rows) is changing .In this paper, we study some effective approaches tothe sign and singularity testing and to the determinantcomputation. We view this classical and modern topicas a good occasion for demonstrating the power ofcombining algebraic and numerical computational techniques.Let us first recall some previous work.%%Motivated by the computational geometry%%applications, one may assume the entries%%of $A$ defined with a high precision, for%%smaller $n$, $n \leq 10$, which is needed%%to ensure the correct answer in the cases%%where det $A$ is close to 0.%%, and very high precision computation is generally%%required already for $n=3$ and $n=4$ if we just%%rely on the straightforward computation of det $A$%%without divisions and represent numbers with%%fixed or floating point. In the case whereIf $n=2$ or $n= 3$, one may try tocompute $\det A$ (so as to output its sign)by using no division and based on its decomposition$\det A = \sum_{j=0}^{n-1} a_{0,j} (-1)^{j-1} A_{0,j}$,where $A_{0,j}$ is the minor of $A=(a_{i,j})$ complementingthe entry $a_{0,j}$. In this case, however, one eitherriscs to lose correct answer due to rounding errors orgenerally needs to double or triple the input precision,which can be prohibitive for some computer implementations.In this and other cases, one may avoid involving higherprecision values by performingelementary transformations of $A$ (that is, byinterchanging its rows and/or columns and by replacinga row or a column by its linear combination with otherrows or columns of $A$). Successful work in this directioncan be traced back to [R52], [E67], and [B68], where thecomputations were performed in rational arithmetic, withno roundoff errors. The recent paper [ABDPY] shows how tobound the precision of computations more effectively inthe case where $n \leq 3$ (the method remains fairlyeffective for $n=4$) and where the objective is not thecomputation of $\det A$ but only of its sign.The algorithm of K.E. Clarkson [C] competes withone of [ABDPY] for $n \leq 4$ and supercedes it forlarger $n$. It relies on the application of the modifiedGram-Schmidt algorithm for the $QR$ (orthogonal)factorization of $A$, which is known to lead to quiteeffective approximate solution and which Clarkson hasextended to the provably correct computation of the sign of thedeterminant. As an example of distinct approaches,we will cite the algorithms that test some classes ofgeometric degeneracies by relying on small (linear)perturbations of the input parameters (cf. [ECS]), andwe refer the reader to [Y], on further information andbibliography.A simple sketchy argument [D] suggests that, under somerestricted model of computing, even the singularitytest (verifying whether det $A = 0$) has the sameasymptotic arithmetic complexity as the evaluation ofdet $A$. Furthermore, under the straight-line programmodel (with no branchings), the evaluation of det $A$and $A^{-1}$ has the same asymptotic arithmeticcomplexity. (The proof of the latter property relieson the partial derivative theorem [L], [BS].)These considerations suggest that decreasing the arithmetic costof computing the sign of $\det A$ is a quite hard problem.\subsection{\bf Our study. Numerical and algebraic approaches.}Our main goal in the present paper is the easier task of decreasing thebit-precision required in algorithms for computing $\det A$or its sign (without blowing up the arithmetic complexitybounds), and such a decrease is highly  importantfor practical applications to computational geometry.We do not confine our study toa single algorithm but consider a broad range ofvarious techniques of numerical linear algebra and rationalalgebraic computations(with no errors).In particular, numerical algorithms for computing variousfactorizations of the input matrix $A$ [that is, itsorthogonal ($QR$) and its triangular ($PLUP_1$, $PLU$and $LU$) factorizations] seem to be the most effectivebasis for computing the sign of $\det A$ providedthat $|\det A|$ is large enough, relatively to computerprecision. Indeed, the sign of $\det A$ can be immediatelyread off from the signs of the diagonal entries ofthe computed triangular matrices, that is, either $R$ or $L$and $U$, respectively.To verify if this output is correct, we estimate the outputerror caused by the accumulation of the rounding errors ofour floating point computation with finite precision.This enables us to obtain  boththe product $d$ of the diagonal entries of $R$ or $L$ and $U$and a positive $e$ such that$\det A$ is bracketted in the interval%%%$[d-e, d+e]$, that is,$d-e \leq \det A \leq  d+e$.Generally, for numerical computations, such goalscan be achieved by means ofinterval analysis (see e.g. [KM]), but in our case,Wilkinson's techniquesof backward error analysis and their  extensionsgive us smaller bounds $e$;moreover, we analyze various algorithms for factorization of$A$, so as to  be able to decrease $e$.If $e<|d|$, then $d$ and$\det A$ have the same sign, and our main problem of thesign computation is solved.It remains to compute the sign of $\det A$ in the casewhere $e \geq |d|$, that is, where\begin{equation}\label{equ1.1}|\det A| < |d| +e \leq  2e.\end{equation}According to the known estimates, based on the backwarderror analysis and recalled in our sections \ref{sec4} and\ref{sec6}, $e$ is proportional to the upper bound $\epsilon = 2^{-\beta}$on the magnitudes of the relative rounding errors, where $\beta$ denotesthe computer precision. Thus, if the inequalities(\ref{equ1.1}) hold, $|\det A|$tends to be substantially smaller than its a priori upperbound $\|A\|^n$. (The numerical algorithms forcomputing $\det A$ or its sign, the analysis of thesealgorithms, and the error estimates are subjects of Part I, sections \ref{sec3}--\ref{sec6}.)Now, if $|\det A|$ is not large, then our problem of obtainingthe sign of $\det A$ by using computations witha lower precision is more easily treatedby means of algebraic methods, which we study in sections\ref{sec7}--\ref{sec15}. In principle, the latter methodscan be applied based on any upper bound on $|\det A|$ but become moreeffective in the cases where $|\det A|$ is smaller,thus complementing the numerical approach. (The algebraic methodsare applied to integer matrices; a real martix canbe turned into an integer matrix via its scalingafter rounding-off its entries.) Indeed, if thecomputed value $d$, approximating$\det A$ numerically, within an absolute error bound $e$, turns out to have a relatively smallmagnitude $|d|$, thenit suffices to compute $(\det A)$ mod $M$, forany integer $M$ exceeding $2(e+|d|)$ (seefact \ref{fac7.1} in section \ref{sec7}).This only requires computations with aprecision of $\lceil {\log_2 M} \rceil$ or fewer bits.Furthermore,  assuming that an upper bound $M$ on$2|\det A| + 1$ is available,we study two customary algebraic approaches to computing$(\det A) \bmod M$, which both are performed with a much lowerprecision than $\lceil \log_2 M\rceil$, except for their finalstages, involving fewerarithmetic operations. Actually, our approach basedon the Chinese remainder theorem and presented insections \ref{sec7}--\ref{sec11} requires lowerprecision computations even at the final stage, providedthat our task is restricted either to testing matrix singularity(see section \ref{sec8}) or,more generally, to determining whether$|\det A|$ is bounded from above by a fixed moderately largenumber and, if so, computing $\det A$ (see section 9).Moreover, even inthe general case of any input matrix $A$, we decrease the precisionof computing at the final stage,at the price of a respectivemoderate increase of the arithmetic cost (see theorem \ref{the11.1}in section \ref{sec11}).This is achieved by means of replacing arithmetic operations with "long"integers by the same operations with associated polynomials having"short" integer coefficients; such a process reverses the knownmethod ofbinary segmentation (see e.g. [BP], pages 277--279, or [P92a]).In section \ref{sec11}, we show how to make all operations withassociated polynomials "linear"; the class of such operations,besides additions andsubtractions, only includes multiplications and divisionsby "short" integers, rather than by polynomials.By avoiding more complicated codes and subroutines for dealingwith nonlinear multiprecision arithmetic, we keepour algorithms for the sign of $\det A$more readilyaccessible for practical implementation.In sections \ref{sec12}--\ref{sec15}, we specify algorithms thatcompute $(\det A) \bmod M$ by using a low precision ofcomputing  and by relying on $p$-adic (Newton-Hensel's) lifting,rather than on the Chinese remainder theorem.For matrices $A$ having larger norms, this approach alsogenerally involves multiplications of pairs of "longer"integers or of their associated nonconstant polynomials,besides linear operations.Among the possible alternative algorithms, we briefly examine(in section \ref{sec16}) the straightforward algorithms amelioratedby means of our approach of section \ref{sec11}; then againonly linear operations with the associated polynomialsare involved in this case. The resulting amelioration may make thealgorithms practically competitive for small $n$.%%%sections \ref{sec7}--\ref{sec14}\subsection{\bf Further comments.}Our techniques of sections 11, 13 and 14 for decreasing theprecision required in the Chinese remainder and $p$-adiccomputations and also our refinement of the estimates of[P85], [P87] for the probability that $(\det A) \bmod p$vanishes for a random prime $p$ provided that $\det A \neq 0$(cf. appendix A) may be of some independent interestfor the designers of algebraic algorithms.The power of the proposed techniques can be accentuatedin their combination with other known ones. In particular,the techniques of [C] can be easily adjusted to thetask of numerical factorization, which we face in sections\ref{sec3}--\ref{sec6}, though we prefer to use Gaussian eliminationwith pivoting, the Householder and/or Givens algorithms tousing the modified Gram-Schmidt algorithm (see our remarks \ref{rem4.2} ofsection \ref{sec4} and \ref{rem6.4} of section \ref{sec6});the techniques of [ABDPY] naturallycomplement ours (in the case of smaller $n$) as an additionalmeans of filtering off the case of absolutely largedeterminants, and a substantial acceleration of numericalfactorization of $A$ and of algebraic computationof $(\det A)$ mod $M$ can be achieved by means of parallelprocessing. [See our remark \ref{rem10.1} in section \ref{sec10}on parallel acceleration of algebraic algorithms formatrix factorizations, see [GL] and [GPS] onparallelization of numerical factorization algorithms,and observe thatat least the first two stages of each of our algorithms\ref{alg7.1}, \ref{alg8.1} and \ref{alg9.1} can be performedconcurrently in $i$, for $i=1, \ldots, k$, by using $k$ processors.]Presenting the estimates for the complexity of matrixfactorizationand of performing our other algorithms, we show theseestimates either exactly or in the form $cn^3 + O(n^2 )$or $cn^3 + O(n^2\log n )$, where we specify the constant $c$.[This is more precise than, for instance, in [C] and [ABDPY],where such bounds have been given in the form $O(n^3)$.]The form $cn^3 + O(n^2 )$ with a specified $c$ is customaryin the field of matrix computations [GL], and we have ignoredthe possible theoretical improvement of these estimates basedon fast matrix multiplication (cf. [BP], [GL]).We have estimated the sequentialarithmetic cost of the presented tests and algorithms, but forseveral of them, effective parallel implementations are wellknown [BP], [GL], [GPS]. For those applications to computationalgeometry where $n$ is small, the asymptotic computationalcomplexity bounds areonly partially informative, however, and the decisivecomparison of the proposed algorithms with each othershould rely on experimental tests.For convenience, we will formally study the evaluation of $\det A$,even though we actually care mostly about computingits sign, due to the major applications to computational geometry.\subsection{\bf Organization of the paper.}After some preliminaries of section \ref{sec2},we present and analyze numerical algorithms for approximating$\det A$ based on numerical factorizations of $A$(see sections \ref{sec3}--\ref{sec6}).In sections \ref{sec7}--\ref{sec9}, we describe the algebraic approachbased on the Chinese remainder theorem,  and  insections \ref{sec12}--\ref{sec15}, we rely on $p$-adic(Newton-Hensel's) lifting.In section 10, we describes the auxiliary computation of$\det A$ modulo a fixed prime,  required in sections\ref{sec7}--\ref{sec9}, \ref{sec12}, and \ref{sec15}.We show our techniques for decreasing theprecision of the Chinese remaindercomputations in section \ref{sec11} and of $p$-adic computations in sections\ref{sec13} and \ref{sec14}.In section \ref{sec16}, we examine the lower precision modificationsof the staightforward algorithms (for smaller $n$).Appendix A contains some estimates for theprobability that reduction modulo a prime chosen at randomin a fixed interval makes thedeterminant of a nonsingular integer matrix vanish. Parts I and IIand appendix A are due to the first author and represent his 1996 revision of his unpublished report of 1994 ``Combining Algebraic, Numericaland Randomization Techniques for the Determinant Computation".In appendix B, the three authors present their numerical experiments forPart I. \ Improved algebraic solution algorithms and the results of the respective tests can be found in [BEPP].{\bf Acknowledgements.} The importance of the problemsstudied in this paper was brought to my attention byKen Clarkson and Franco Preparata. Ken Clarkson andJohn Hobby also described the customary requirementsto the input matrices encountered in the applicationsto computational geometry, and Ken Clarkson made severaluseful comments to the first draft of this paper.Bas Braams and Richard Pollack gave me some usefulinformation about the computational geometry applicationsin the case of matrices of larger sizes. The discussionswith Herv\'{e} Br\"onnimann and, more briefly, withOlivier Devillers on implementation and applicationsof determinant sign algorithms and Herv\'{e}'s comments on theworks [C] and [ABDPY] were very pleasant and informative.%%%%Helpful criticism by the referees was also a%%%%springboard for my additional work on this paper.\vspace{.1in}\noindent{\bf PART I. NUMERICAL COMPUTATION APPROACH.}\vspace{.1in}\noindent\section{\bf Basic matrix factorizations.} \label{sec2}\setcounter{equation}{0}$\det A$ and its sign for an $n \times n$ integer,rational or real input matrix $A$ can be immediately obtained from the factorizations:\begin{equation}\label{equ2.1}A=LU\  ,\end{equation}\begin{equation}\label{equ2.2}A=PLU\  ,\end{equation}\begin{equation}\label{equ2.3}A=PLUP_1\  ,\end{equation}\begin{equation}\label{equ2.4}A=QR\  ,\end{equation}where $Q$ is a real $n \times n$ orthogonal (unitary) matrix,such that\begin{equation}\label{equ2.5}Q^TQ = I\  ,\end{equation}$L^T$, $U$ and $R$ are $n \times n$ upper triangularmatrices, and $P$ and $P_1$ are permutation matrices,whose application to any vector $v$amounts to some fixed permutations of its coordinates.Here and hereafter, $W^T$ denotes the transpose of amatrix or a vector $W$, and $I$ and $O$ denotethe identity matrices, $Iv=v$ for any vector $v$, and thenull matrices of appropriate sizes. The computation of $\det A$ based on(\ref{equ2.2})--(\ref{equ2.5})relies on the following well-known facts:\begin{fact}\label{fac2.1}$\det (SV) = (\det S)\det V$ for any pair of$n \times n$ matrices $S$ and $V$.\end{fact}\begin{fact}\label{fac2.2}$\det T =  \displaystyle{ \prod_{i=0}^{n-1} } t_{i,i}$for any $n \times n$ triangular matrix$T=( t_{i,j}, \ i, j=0, \ldots, n-1)$.\end{fact}\begin{fact}\label{fac2.3}$\det P=1$ or $\det P=-1$ for any permutation matrix $P$;$\det I = 1$.\end{fact}\begin{fact}\label{fac2.4}$\det W=\det W^T$ for any square matrix $W$.\end{fact}\begin{corollary}\label{cor2.5}$\det Q=1$ or $\det Q =-1$ for any real square matrix $Q$satisfying (\ref{equ2.5}).\end{corollary}%%\begin{remark}\label{rem2.1}%%The case of complex matrices $A$ can be treated similary,%%by using complex orthogonal (unitary) matrices $Q$%%satisfying the matrix equation%%$Q^UQ=I$, where $Q^T$ denotes the Hermitian transpose of $Q$.%%\end{remark}%%\section{\bf Solution based on Gaussianelimination.}\label{sec3}\setcounter{equation}{0}Based on the equation (\ref{equ2.3}) andfacts \ref{fac2.1}--\ref{fac2.3}, wedevise the following algorithm for computing $\det A$ and/orits sign:\begin{algorithm}\label{alg3.1} $ \   ${\bf Input:} an $n \times n$ matrix $A$.{\bf Output:} $\det A$.{\bf Computations.}\begin{enumerate}\itemCompute the matrices $P, \  L, \ U$, and $P_1$ satisfying(\ref{equ2.3}).\itemCompute $\det P, \  \det L, \ \det U$, and $\det P_1$.\itemCompute and output \\ $\det A = (\det P)( \det L)(\det U)(\det P_1)$.\end{enumerate}\end{algorithm}Correctness of the algorithm immediately follows from facts\ref{fac2.1}--\ref{fac2.3}.If we only need to output thesign of $\det A$, then we only need to compute the signsof $\det P, \  \det L, \ \det U$, and $\det P_1$ at stage 2and to multiply these signs at stage 3.The computation of the matrices $P, \  L, \  U$, and $P_1$ atstage 1 goes by means of routine application of {\em Gaussianelimination with complete pivoting}, which,as a by-product, records the permutations defined by the matrices$P$ and $P_1$.The values ($1$ or $-1$)  of $\det P$ and $\det P_1$ can beimmediately  recovered from these records. Furthermore, thematrix $L$ is a unit lower triangular  matrix, whosediagonal is filled with ones, so that $\det L = 1$.The computation of $\det U$ only requires $n-1$ multiplications,due to fact \ref{fac2.2}. The overall computational cost ofperforming algorithm \ref{alg3.1}is dominated  by the cost of performing its stage 1:$\sum_{i=1}^{n-1} i^2=(2n-1)(n-1)n/6$ multiplications,as many additions/subtractions, as many comparisons, and$\sum_{i=1}^{n-1} i = (n-1)n/2$ divisions.{\em Gaussian elimination with partial pivoting} computes the factorization(\ref{equ2.2}) [that is, (\ref{equ2.3}) for $P_1=I$] and uses$(n-1)n/2$ comparisons, rather than $(2n-1)(n-1)n/6$ (compareremark \ref{rem4.2} in the next section).\section{\bf Numerical implementation of the solutionbased on Gaussian elimination.}\label{sec4}\setcounter{equation}{0}In its numerical implementation with a finite precision,algorithm \ref{alg3.1} computes the triangular  matrices$\tilde{L}=L+E_L$ and $\tilde{U}=U+E_U$, as well as the permutationmatrices $\tilde{P}$and $\tilde{P_1}$, where $E_L$ and $E_U$ denotethe matrices of the pertubations of$L$ and $U$, due to accumulation of rounding errors.(Since rounding errors are usuallysmall, we usually have $\tilde{P}=P$and $\tilde{P_1}=P_1$, but we do not need to assume this.)\begin{definition}\label{def4.1}Write\begin{equation}\label{equ4.1}A'=\tilde{P}^{-1}A\tilde{P_1}^{-1}, \\tilde{A}'= A'+E_{A'}=\tilde{L}\tilde{U}\end{equation}and  let $e'$, $a$, $\tilde{\ell}$ and$\tilde{u}$ denote the maximumabsolute values of the entries of the matrices $E_{A'}$,$A, \ \tilde{L}$ and$\tilde{U}$, respectively.\end{definition}Our next goal is to estimate $e'$ from above, assumingfloating point binary arithmetic with rounding to $\beta$ bits. This means  the upper bound\begin{equation}\label{equ4.2}\epsilon = 2^{-\beta}\end{equation}on the magnitudes of the relative rounding errors (sothat $\epsilon$ is an upperbound on the {\em unit roundoff}, also called the {\em machineepsilon}).\begin{theorem}\label{the4.1} {\rm (cf. e.g.  [CdB], page 181,and [SB],pages 191--196, for derivation and minor refinements). }For a matrix $W=(w_{i,j})$, let $|W|$ denote the matrix$(|w_{i,j}|)$ obtainedby replacing each entry $w_{i,j}$  with its absolute value$|w_{i,j}|$. Then, under(\ref{equ4.1}), we have$$ |E_{A'}| \leq (|A'|+|\tilde{L}| |\tilde{U}| )\epsilon. $$\end{theorem}\begin{corollary}\label{cor4.2}\begin{equation}\label{equ4.3}e'\leq e^+=(a + n\tilde{\ell}\tilde{u})\epsilon .\end{equation}\end{corollary}In this paper, we will also express some estimates in termsof matrix norms,$\|W\|_q$, $q=1,2, \infty$, satifying the following basicrelations, where$W=(w_{i,j})$:\begin{equation}\label{equ4.4}\|W\|_1 = \max_j \sum_i |w_{i,j}|, \  \|W\|_{\infty} =\max_i \sum_j |w_{i,j}|,\end{equation}\begin{equation}\label{equ4.5}(1/n)\|W\|_q  \leq  \max_{i,j} |w_{i,j}| \leq \|W\|_q, \ q=1,2,\infty,\end{equation}\begin{equation}\label{equ4.6}\|W\|_2^2  \leq  \|W\|_1\|W\|_{\infty}, \   \|W\|_2 \leq \sqrt n \|W\|_q, \ q=1, \infty\end{equation}(see [GL], pages 56--58).The norm $\|W\|_2$ is actuallyneeded only in section \ref{sec6}.We obtain from (\ref{equ4.5}) that$e' \leq \|E'_A\|_{\infty}$ and recall the followingresult, in whose statement we use an upper bound$\epsilon$ of (\ref{equ4.2}):\begin{theorem}\label{the4.3}{\rm ([FM], chapter 21). }\\\begin{equation}\label{equ4.7}e' \leq e_+=\|E'_A\|_{\infty} \leq n^2 a^+\epsilon,\end{equation}where $a^+$ denotes the maximum absolute value of allthe entries of thematrices $A^{(i)}$ overwriting $A'$ and computed in theprocess of Gaussianelimination, which reduces $A'$ to the upper triangular form.\end{theorem}\begin{remark}\label{rem4.1} Ignoring the termsproportional to $\epsilon^k, \ k\geq 2$, we mayimprove (\ref{equ4.7})as follows: $e' \leq 3(n-1)\epsilon a^+ + O(\epsilon^2)$(see [DB], page 178). The latter bound and both (\ref{equ4.3}) and (\ref{equ4.7}) area posteriori  bounds:  besides the input data, they depend on the values$\tilde{\ell}, \ \tilde{u}$ or $a^+$, which are available as by-productsof our computations.In the context of our task, a priori bounds,depending only on the input data may only be needed for theoreticalanalysis. Such bounds can beobtained based on Wilkinson's estimate (cf. [Wil]):$$a^+ \leq an^{1/2} \prod_{k=1}^{n-1} (k+1)^{1/(2k)} <1.8 n^{({\rm \ell n} \  n )/4} a, \ a=\max_{i,j} |a_{i,j}|.$$Combining this estimate with theorem \ref{the4.3} implies that$$e' \leq  e_+ =\epsilon an^{2.5}\prod_{k=1}^{n-1} (k+1)^{1/(2k)} <1.8 \epsilon an^{2 + ({\rm \ell n}\  n)/4}, \ a=\max_{i,j} |a_{i,j}|.$$The above a priori bounds, however, are overly pessimistic accordingto the extensive practical experience with Gaussianelimination, and this has even motivated the conjecture that, actually,for all matrices $A$, we have\begin{equation}\label{equ4.8}a^+=O(an), \  \  e'=O(\epsilon a n^3),\end{equation}in the case of complete pivoting.\end{remark}Now, we wish to estimate the magnitude of the perturbationof $\det A$ and$\det A'$ caused by the perturbation of $A'$. To simplifythe notation, wewill assume that $P=P_1=I$, $A'=A$ and will write\begin{equation}\label{equ4.9}e=e', \  \   e_d = \det (A+E_A)-\det A\end{equation}(cf. definition \ref{def4.1}).We will need the following simple and/or well-known estimates:\begin{fact}\label{fac4.4}$|\det A| \leq \|A\|_q^n, \  q=1,2, \infty$.\end{fact}\begin{fact}\label{fac4.5}$\|B\|_q \leq \|W\|_q, \  q=1,\infty$,  if$|B| \leq  |W|$; \\$\|B\|_q \leq \|W\|_q$, $q=1,2,\infty$,  if $B$is a submatrix of $W$.\end{fact}\begin{fact}\label{fac4.6}$\|A+E_A\|_q \leq \|A\|_q + n e$, for  $q=1, 2, \infty$;$\|A+E_A\|_2^2 \leq \|A\|_2^2 + n e$.\end{fact}By using these results, we estimate $e_d$ of (\ref{equ4.9}).\begin{proposition}\label{pro4.7}$|e_d| \leq  (\|A\|_q + ne)^{n-1} n^2 e $, for $q=1, 2, \infty$.\end{proposition}{\bf Proof.}  Represent $\det(A+E_A)$ as the sum ofmonomials$\pm \tilde{a}_{i_0,0} \tilde{a}_{i_1,1}\ldots \tilde{a}_{i_{n-1},n-1}$,where $\tilde{a}_{i_g,g} =a_{i_g,g} +e_{A,i_g,g} ,\  g=0,1,\ldots, n-1$.Rewrite such a monomial as a telescopic sum,$$\pm [ (\tilde{a}_{i_0,0} - a_{i_0,0}) \tilde{a}_{i_1,1}\ldots \tilde{a}_{i_{n-1},n-1} +  \\a_{i_0,0}(\tilde{a}_{i_1,1} - a_{i_1,1}) \tilde{a}_{i_2,2}\ldots \tilde{a}_{i_{n-1},n-1} + \ldots +$$$$ a_{i_0,0}a_{i_1,1} \ldots a_{i_{n-2},n-2}(\tilde{a}_{i_{n-1},n-1} - a_{i_{n-1},n-1}) +a_{i_0,0} a_{i_1,1} \ldots a_{i_{n-1},n-1}].$$Summation of the respective terms of this sum over allthe subscripts of the monomials of $\det (A+E_A))$gives us the following equation:$$\det (A+E_A) = \sum_{i,j} (\tilde{a}_{i,j} - a_{i,j})\tilde{D}_{i,j} +\det A,$$where $\tilde{D}_{i,j}$ denotes the determinantof the $(n-1) \times (n-1)$ matrix obtained byreplacing the first $j-1$ columns of $A$ with thefirst $j-1$ columns of $\tilde{A} = A+E_A$ and by deletingthe $i$-th row and the $j$-th column of the resulting matrix.Therefore,$$|e_d|=|\det (A+E_A)-\det A| =|\sum_{i,j} (\tilde{a}_{i,j} - a_{i,j}) \tilde{D}_{i,j}| \leq\sum_{i,j} |e_{A,i,j}| |\tilde{D}_{i,j}|$$$$\    \  \  \leq  \   e \sum_{i,j} |\tilde{D}_{i,j}|\leq n^2 e \max_{i,j}|\tilde{D}_{i,j}|.$$ By combining facts \ref{fac4.4}--\ref{fac4.6}, we obtain that$|\tilde{D}_{i,j}| \leq (\|A\|_q + ne)^{n-1}$ for $q=1, 2, \infty$.Substitute the latter estimates into theabove bound on $|e_d|$and yield proposition \ref{pro4.7}.  \hfill $\Box$Combining proposition \ref{pro4.7} with the bound (\ref{equ4.7})enables us to extend algorithm \ref{alg3.1} as follows (cf. remark \ref{rem4.2} of this section):\begin{algorithm}\label{alg4.1} $ \  \  ${\bf Input:} an $n \times n$ real matrix $A$ and apositive $\epsilon$.{\bf Output:} a pair of real numbers $d_-$ and $d_+$ such that$d_- \leq \det A \leq d_+$.{\bf Computations.}\begin{enumerate}\itemApply algorithm \ref{alg3.1} in floatingpoint arithmetic with unit roundoff bounded by $\epsilon$.Let $\tilde{U}=U+E_U$ denotethe computed upper triangular matrix, approximating the factor$U$ of $A$ [from (\ref{equ2.2}) or (\ref{equ2.3})].\itemCompute the upper bound $e_+$ of (\ref{equ4.7}) on $e=e'$of  definition \ref{def4.1}.\itemSubstitute $e_+$ for $e$ and $\min ( \|A\|_1, \  \|A\|_{\infty}, \(\|A\|_1 \|A\|_{\infty})^{1/2})$ for $ \|A\|_q$in proposition \ref{pro4.7} and obtain an upper bound $e_d^+$ on $|e_d|$.\itemOutput the values $\det \tilde{U} - e_d^+$and $\det \tilde{U} + e_d^+$.\end{enumerate}\end{algorithm}Correctness of the algorithm follows from theorem \ref{the4.3}and proposition \ref{pro4.7}.Due to (\ref{equ4.3}), algorithm \ref{alg4.1} can be modifiedwith $e^+$ replacing $e_+$ at its stages 2 and 3. In bothversions, the computations at stages 2--4 involve $O(n^2)$arithmetic operations (hereafter referred to as {\em ops})and comparisons [cf. (\ref{equ4.3}), (\ref{equ4.4}), (\ref{equ4.7})],which is dominated by thearithmetic cost of performing stage 1.We may obtain the sign of $\det A$ as the sign of$\det \tilde{U}$ provided that\begin{equation}\label{equ4.10}e_d^+ <|\det \tilde{U}|.\end{equation}\noindentOtherwise, some additional computations are needed. They can besimplified  (see section \ref{sec7}) by using the bound$$|\det A| \leq e_d^+ + |\det \tilde{U}|,$$which turns into\begin{equation}\label{equ4.11}|\det A| \leq 2 e_d^+\end{equation}unless (\ref{equ4.10}) holds [compare (\ref{equ1.1})].For smaller $\epsilon$, the bound (\ref{equ4.11}) issubstantially stronger than the bounds of fact \ref{fac4.4}(compare proposition \ref{pro4.7}).\begin{remark}\label{rem4.2} Algorithm \ref{alg4.1}performed with partial (rather than complete) pivotingnumerically computes the factorization(\ref{equ2.2}).The advantage is in saving$\sum_{i=1}^{n-1}(i^2-i) = (n-2)(n-1)n/3$ comparisons.The disavantage is a plausible substantial growth of the upperbounds on $e=e'$ and $e_d$. Namely, Wilkinson's a prioriupper bound on $a^+$ of remark \ref{rem4.1} does not holdanymore and is replaced by a much weaker a priori upperbound, $a^+ \leq 2^{n-1} a, \ a=\max_{i,j}|a_{i,j}|$.The latter bound has been reached for some specially concoctedmatrices $A$, but extensive experience of performingGaussian elimination in worldwide computational practicesuggests  that, as a rule, the value $a^+$ remains boundedfrom above by the value $cna$, for some moderate constant $c$, evenin the variant with partial pivoting [GL]. Thissuggests that algorithm \ref{alg4.1} should be actually appliedwith using partial pivoting.In the unlikely case of excessively large $a^+$, we stillmay repeat the computations withcomplete pivoting or we may shift toalgorithms \ref{alg5.1} and/or\ref{alg6.1} of the next sections, which compute$QR$ factorization of $A$.For the latter algorithms, the output error bounds that are veryclose to (\ref{equ4.8}) have been formally provenand not only conjectured  [see (\ref{equ6.4}) insection \ref{sec6}]. Furthermore, a modificationof these algorithms (based on using Givens rotations)is highly effective for computing (the sign of)$\det A$ for a dynamically updated matrix $A$, which isunmatched if we rely on algorithms \ref{alg3.1} and\ref{alg4.1} (see the end of the next section).Those readers for who these advantages, however,  are notimportant may skip the next two sections andgo directly to section \ref{sec7}.\end{remark}\section{\bf Solutions based on $QR$ factorization.}\label{sec5}\setcounter{equation}{0}Next, based on the equations (\ref{equ2.4}) and (\ref{equ2.5}),fact \ref{fac2.1}, \ref{fac2.2}, and corollary \ref{cor2.5},we devise the following algorithm, as an alternative toalgorithm \ref{alg3.1}:\begin{algorithm}\label{alg5.1} $ \  \   ${\bf Input:} an $n \times n$ matrix $A$.{\bf Output:} $\det A$.{\bf Computations.}\begin{enumerate}\itemCompute a matrix $Q$ satisfying (\ref{equ2.4}) and(\ref{equ2.5}) and represented either by its entries or bya sequence of its orthogonal factors.\itemCompute $\det Q$.\itemCompute the matrix $R=Q^TA=(r_{i,j})$.\itemCompute $\det A = (\det Q) \prod_i r_{i,i}$.\end{enumerate}\end{algorithm}Correctness of the algorithm immediately followsfrom facts \ref{fac2.1} and \ref{fac2.2}.We will consider two effective customary algorithmsfor computing the $QR$ factorization at stage 1, that is,the Householder and the Givens algorithms (see [Wil],[GL] and compare our remark \ref{rem6.4} in thenext section). The Givens algorithm applied at stage 1always produces a matrix $Q$ with $\det Q =1$ since itcomputes $Q$ as the product of the {\em matrices of Givensrotations}, each having determinant 1 (see the matrices$G$ in the beginning of the next section or see [GL],page 201). Therefore, stage 2 is a by-product of stage 1 inthis case. Let us deduce a similar property in the casewhere we performstage 1 by means of the Householder algorithm, which computes$Q$ as a product of $n-1$ matrices of Householdertransformations.\begin{fact}\label{fac5.1}For a matrix $Q$ of (\ref{equ2.4}), (\ref{equ2.5}) computedby means of the Householder algorithm, we have\end{fact}\begin{equation}\label{equ5.1}\det Q = (-1)^{n-1}.\end{equation} %%$${\bf Proof:} Recall that the latter algorithmdefines $Q$ as the product of $n-1$ Householder matrices,$Q=\prod_{i=1}^{n-1}H_i \ , \ \ H_i = I -2v_iv_i^T / v_i^Tv_i$for some vectors $v_i, \ \ i= 1, \ldots, n-1$. Now,(\ref{equ5.1}) follows since\begin{equation}\label{equ5.2}\det H_i = -1 \ , \ \ {\rm for \ all } \ \ i \ .\end{equation}To prove (\ref{equ5.2}), recall that, $H_i^TH_i = I$,so that $(\det H_i)^2=1$. On the other hand,$D(0)=\det (I-2 e_0e_0^T) = -1$, where$e_0^T=(1, 0, \ldots, 0)$, $e_0^Te_0 = 1$.Consider the homotopic transformation$v(t) = e_0 + t(v_i - e_0)$, for $t$ ranging from 0 to 1and observe that $D(t) = \det (I-v(t)v^T(t)/( v^T(t)v(t)))$is a function continuous in $t$, $D^2(t) = 1$ for all $t$,and $D(0)= -1$. Therefore,$D(1) = \det I -v_iv_i^T/v_i^Tv_i = D(0)=-1$,and (\ref{equ5.2}) follows.\hfill $\Box$Stage 3 of algorithm \ref{alg5.1} uses $n-1$ multiplications.Stage 1 is, by far, the hardest. It involves $(4/3)n^3 + O(n^2)$multiplications, about as many additions/subtractions, and$O(n)$ evaluations of the square roots of positive numbersif the Householder algorithm is applied. With the Givensalgorithm, these bounds turn into $2n^3 + O(n^2)$,$n^3 +O(n^2)$, and $O(n^2)$, respectively.This shows some minor advantageof the Householder algorithm, and similarly does thecomparison of the associated bounds on the output perturbationcaused by rounding errors [compare (\ref{equ6.4})and (\ref{equ6.5}) in the next section].Using Givens rotations, however, is highly effective inthe cases where the matrix $A$ is dynamically updated,that is, where some of its rows and/or columns aredynamically replaced by new ones. Givens rotationsenable us to update the $QR$ factorization of $A$ successively, byusing $O(k n^2)$ ops and square root evaluations,  if$O(k)$ rows and/or columns of $A$ are successively replaced by newones (see [GL], sections 12.6.1-- 12.6.3,pages 593--597).This is a special feature of Givens rotations, not sharedby other $QR$ factorization algorithms. Furthermore,dynamic updatings of the $PLU$ factorizationof (\ref{equ2.2})  and of the $PLUP_1$ factorization of(\ref{equ2.3}) are more complicatedand more involved, due to the problem of pivoting [GGMS].On the other hand, the efficacy of the fast updating of the $QR$factorization by means of Givens rotations is limited because,after $s$ successive fast updatings, the implicitly computedmatrix $Q$ is represented as a product of $r=s+n(n-1)/2$Givens matrices, which means recursive accumulation of many roundingerrors [see the bound (\ref{equ6.3}) in the next sectionfor larger $r$].Therefore, one needs to recompute the $QR$ factorizationof the current input matrix from the scratch, after sufficientlymany updatings of $A$, when $s$ and $r$ grow too large.%and (\ref{equ2.5}) lar\section{\bf Numerical computation of the orthogonal factor $Q$(in a factorized form).}\label{sec6}\setcounter{equation}{0}It is convenient to keep $Q$ in the factorizedform, as a product of the matrices of Householdertransformations or Givens rotations [GL], [Wil].Keeping also the Householder matrices in theform $H=I-2vv^T$, for the associated vectors $v$,we preserve the equation $\det Q=(-1)^{n-1}$, evenin the case of numericalimplementation of stage 1, with a finite precision and roundingerrors.Let us examine the perturbation of $Q$ caused by rounding errorsin the case of finite precision implementation of the Givensalgorithm. The matrices of Givens rotations are of the form$G=I+R_{i,k}(c,s)$, where the matrix $R_{i,k}(c,s)$ hasexactly four nonzerosentries,  that  is,  $R_{i,i}(c,s)=R_{k,k}(c,s) = c-1, \R_{i,k}=-R_{k,i}=s$,for two fixed integers $i$ and $k$  and for a pair of nonzero real$c$ and $s$,satisfying $c^2+s^2=1$, so that $G^TG=I$.In a finite precision implementation, $c$ and $s$ appearas $c+e_c$and $s+e_s$, respectively, where $e_c$ and $e_s$ denote therounding errors,\begin{equation}\label{equ6.1}|e_c|< |c| \epsilon, \  |e_s|< |s| \epsilon,\end{equation} and $\epsilon = 2^{-\beta}$ bounds the magnitudes of relativerounding errors [cf. (\ref{equ4.2})].Then, $\det (G^TG) = (\det G)^2 = (c+e_c)^2 +(s+e_s)^2 =(c^2 +s^2)(1+\Delta)=1+\Delta$, where, due to the bounds (\ref{equ6.1}),we have\begin{equation}\label{equ6.2}|\Delta|<(1+\epsilon)^2 -1= 2\epsilon +\epsilon^2.\end{equation}Thus, for small $\epsilon$, $|\Delta|$ is small, and$\det G$ is close to 1.  We will assume that, at  stage 1 ofalgorithm \ref{alg3.1}, we represent$Q$ as the product of $r$ Givensmatrices, for $r\leq (n-1)n/2$, and use a finite precision with the bound$\epsilon$ on the relative rounding errors. Then, thematrix $Q$ is actually replaced by the matrix $Q+E_Q$,  $E_Q$denoting the perturbation caused by the rounding errors,and due to (\ref{equ6.2}),  we have\begin{equation}\label{equ6.3}(1-2\epsilon -\epsilon^2)^r \leq \det(Q+E_Q)\leq (1+2\epsilon +\epsilon^2)^r,  r\leq(n-1)n/2.\end{equation}Let us next estimate the perturbation ofthe matrix $R$ caused by rounding errors.According to Wilkinson's {\em the backwarderror analysis}, the rounding errors of floating point computationcause the same perturbation $E_R$ ofthe matrix $R$ as would have followed if the input matrix$A$ had been perturbed by a certain matrix $E=E_A$ and ifthe subsequent computation of $Q$ and $R$ had beenperformed exactly, with no errors [Wil]. Furthermore,Wilkinson has proved the following upper bounds on $\|E\|_2$,assuming an upper bound $\epsilon$ on the unit roundoff:\begin{equation}\label{equ6.4}\|E\|_2 \leq h \|A\|_2 , \  h < h^+=3.35(n-1)(1+9.01\epsilon)^{n-2}\epsilon,\end{equation}\begin{equation}\label{equ6.5}\|E\|_2 \leq g \|A\|_2 , \  g < g^+=3n^{5/2}(\frac{1+6\epsilon}{1-4.3\epsilon})^{2n-4}\epsilon;\end{equation}the bound (\ref{equ6.4}) or (\ref{equ6.5}) holdwhere the Householder or, respectively, Givensalgorithms have been applied at stage 1 of algorithm \ref{alg5.1}(cf. [Wil], pages 236 and 240, respectively).\begin{remark}\label{rem6.1} The reader is referred to [S77] onfurther refinements of theerror bounds for finite precisioncomputation of the $QR$ factorization of a matrix.\end{remark}\begin{remark}\label{rem6.2} The  bound (\ref{equ6.5})incorporates the effect of the perturbation of the matrix$Q$ caused by rounding andestimated in (\ref{equ6.3}).\end{remark}By combining the bounds (\ref{equ4.5}), (\ref{equ6.4}) and(\ref{equ6.5}) with proposition \ref{pro4.7}, we obtainthe following result:\begin{proposition}\label{pro6.1}Let algorithm \ref{alg5.1} be performed by using thefloating point arithmetic with relative rounding errors (unit roundoff)within $\epsilon$.Let at stage 1 of this algorithm, the Householderor Givens algorithms compute numerically the Householderor Givens factors of $Q$, respectively.Then algorithm \ref{alg5.1} outputs the values$d=(-1)^{n-1}\det \tilde{R}$ or $d=\det \tilde{R}$,respectively, satisfying the following relations:$$ d=\det A + e_d, $$\begin{equation}\label{equ6.6}|e_d | \leq n^2 s \|A\|_2 ( \|A\|_q + ns\|A\|_2)^{n-1}, \q=1, 2, \infty,\end{equation}where $s$ stands either for $h$ of (\ref{equ6.4}) orfor $g$ of (\ref{equ6.5}), respectively.\end{proposition}Now, we may specify a floating point implementation ofalgorithm \ref{alg5.1} as follows:\begin{algorithm}\label{alg6.1} $ \   \  ${\bf Input} and {\bf Output}, as in algorithm \ref{alg4.1}.{\bf Computations.}\begin{enumerate}\itemApply algorithm \ref{alg5.1} by using the Householderalgorithm at its stage 1 and by performing all thecomputations in floating point arithmetic withunit roundoff bounded by $\epsilon$. Let $\tilde{R}=R+E_R$ denotethe computed upper triangular matrix approximatingthe factor $R$ of $A$ in (\ref{equ2.4}).\itemSet $s=h$ and specify the two upper bounds of(\ref{equ6.6}), on $|e_d|$, for $q=1$ and $q=\infty$,by substituting for $\|A\|_2$ its upper bound$\min \{\|A\|_1 \sqrt n , \  \|A\|_{\infty}\sqrt n, \(\|A\|_1 \|A\|_{\infty})^{1/2}\}$,[compare (\ref{equ4.6})]. Of the two computed upperbounds on $|e_d|$,  choose the smaller one anddenote it $e_d^+$.\itemCompute and output the values $\det \tilde{R}-e_d^+$ and$\det \tilde{R} + e^+_d$.\end{enumerate}\end{algorithm}Correctness of the algorithm follows from (\ref{equ4.6})and propositions \ref{pro4.7} and \ref{pro6.1}.By replacing the Householder algorithm bythe Givens algorithm and by substituting$s=g$ in to (\ref{equ6.6}), we arrive at analternative extension ofalgorithm \ref{alg5.1}. In both cases, thecomputation at stages 2 and 3of these extensions only requires $O(n^2)$ops and comparisons [see(\ref{equ4.4})].\begin{remark}\label{rem6.3}Generally,  computing the entries of the matrix$Q$ explicitly gives us no advantage in the context of ourtopic of computing $\det A$.Indeed, such a computation of $Q$ roughly doubles thecomputational cost of performing algorithms \ref{alg5.1}and \ref{alg6.1} but does not affect the computed matrix $\tilde{R}$.\end{remark}\begin{remark}\label{rem6.4}At stage 1 of algorithm \ref{alg5.1}, we could havealternatively appliedthe modified Gram-Schmidt (MGS) algorithm, whichis one of the customary tools for computing the$QR$ factorization. The MGS algorithmrequires $4n^3+O(n^2)$ ops in order to compute thefactor $R$ and also needs some additionalcomputation in order to find out if $\det Q$equals $1$ or $-1$.In the case of implementation of the MGS algorithm infloating point arithmetic, the extra factor of ${\rm cond}_2A$(versus the Householder and Givens algorithms) appears in the upperestimate for the 2-norm of the perturbation of the matrix $Q$caused by rounding to a fixed finite precision(see [GL], page 219, or [B67] and note that ${\rm cond}_2 A$ isvery large for nearly singular matrices $A$),which is the main deficiency of the MGS algorithm.\end{remark}\vspace{.1in}\noindent{\bf PART II. ALGEBRAIC COMPUTATION APPROACH.}\vspace{.1in}\noindent\section{\bf Recovery of $\det A$ from $\det A$  modulo an integer$M$ based on the Chinese remainder theorem.}\label{sec7}\setcounter{equation}{0}Hereafter, $a \bmod M$ (for two integers $a$ and $M>1$)denotes an integer$b$ satisfying the relations $0 \leq b \leq M$, $M$ divides $b-a$.In this and the next sections, we will rely on the followingobservation:\begin{fact}\label{fac7.1}For  any pair of integers $d$ and $M>2|d|$, we have\begin{equation}\label{equ7.1}d= (d \bmod M)  \  \  \   if \   d \bmod M < M/2, \  \d =-M + (d \bmod M)  \    \  \   otherwise.  \\\end{equation}\end{fact}Suppose that the input matrix $A$ is filled with  integers.If we know an upper bound $d^+ \geq |\det A|$,we may choose  any integer $M$ exceeding  $2d^+$ andthen compute $(\det A) \bmod M$and recover $\det A$ by applying (\ref{equ7.1}).In section 10, we will review the computation of $(\det A)\bmod M$ for a fixed integer $M>1$. Computations modulo $M$by these algorithms can be performed with theprecision of $  \lceil  \log_2 M\rceil $ bits, butwe wish to perform mostof them with a lower precision. In this and the two nextsections, we will rely on the Chinese remainder theorem [AHU],  [BM].\begin{theorem}\label{the7.1} {\rm (Chinese remainder theorem).}Let $m_1, \ldots, m_k$ denote $k$ pairwise relativelyprime integers (say, $k$ distinct primes),\begin{equation}\label{equ7.2}m_1 > m_2 >  \ldots >m_k>1.\end{equation}Let $D$ denote an integer satifying\begin{equation}\label{equ7.3}0\leq D < M =m_1  m_2   \ldots m_k.\end{equation}Let\begin{equation}\label{equ7.4}r_i = D \bmod  m_i,\end{equation}\begin{equation}\label{equ7.5}M_i= M/m_i,  \  v_i=M_i \bmod m_i, \w_i v_i=1 \bmod m_i, \  i=1, \ldots, k.\end{equation}Then $D$ is a unique integer satisfying (\ref{equ7.3}) and(\ref{equ7.4}); furthermore,\begin{equation}\label{equ7.6}D=(\sum_{i=1}^k M_ir_iw_i)\bmod M.\end{equation}\end{theorem}\begin{algorithm}\label{alg7.1},  computation of $(\det A) \bmod M$based on the Chinese remainder theorem.{\bf Input:} an integer matrix $A$, a black box algorithmthat computes $ (\det A) \bmod m$ for any fixedinteger $m>1$, $k$ integers$m_1, \ldots, m_k$ satisfying  (\ref{equ7.2}) and pairwise relatively prime,and $M=m_1 m_2\ldots m_k$.{\bf Output:} $(\det A) \bmod M$.{\bf Computations.}\begin{enumerate}\itemCompute $r_i =( \det A) \bmod m_i, \  i=1, \ldots, k$.\itemCompute the integers $M_i, \  v_i$ and $w_k$ of(\ref{equ7.5}), for $i=1, \ldots, k$.\itemCompute and output $D$ of (\ref{equ7.6}).\end{enumerate}\end{algorithm}Correctness of the algorithm immediately follows fromtheorem \ref{the7.1}.We compute the values $M_i$, $v_i$ by using (\ref{equ7.5}).The values $w_i$ are obtained by means of applying theEuclidean algorithm to the pair of $v_i$ and $m_i$, $i=1, \ldots, k$(cf. [AHU], [BM]).We need $O(kn^3)$ arithmetic operations (ops) at stage 1,$O((k\log m_1) \log \log m_1)$ ops at stage 2, and $O(k \log^2 k)$ opsat stage 3.The computations are performed with the precision of at most$\lceil  \log_2 m_1 \rceil$ bits at stage 1 and of at most$\lceil  \log_2 M \rceil $ bits at stages 2 and 3.Furthermore, in section 11 we will show how to decrease the latterprecision (at stages 2 and 3) by using a little more ops.\section{\bf Matrix singularity test.}\label{sec8}\setcounter{equation}{0}Next, we will show a modification of algorithm \ref{alg7.1},  for testingmatrix singularity, where stages 2 and 3 are replaced by much simplercomputations.\begin{algorithm}\label{alg8.1},  matrix singularity test.{\bf Input,} as in algorithm \ref{alg7.1}, where\begin{equation}\label{equ8.1}M > |\det A|.\end{equation}{\bf Output:} SINGULAR if $\det A= 0 \bmod M$, NONSINGULAR otherwise.{\bf Computations:}\begin{enumerate}\itemas in algorithm \ref{alg7.1}.\itemOutput SINGULAR if $r_i=0$ for $i=1, \ldots, k$;output NONSINGULAR  otherwise.\end{enumerate}\end{algorithm}If $r_i=(\det A) \bmod m_i \neq 0$ for some $i$, then, clearly,$\det A \neq 0$, and the output of the algorithm is correct.If $r_i=0\bmod m_i$ for all $i$, then $D=0$  due to(\ref{equ7.6}), where $D=(\det A) \bmod M$,by the virtue of  theorem \ref{the7.1}.Due to (\ref{equ8.1}), the latter equations imply that$\det A=0$, and this completescorrectness proof for algorithm  \ref{alg8.1}.  \hfill  $\Box$\begin{remark}\label{rem8.1} Suppose that we have only computed$r_1=(\det A) \bmod m_1$. If $r_1 \neq 0$, then $\det A \neq 0$,so that we arrive at an one-sided singularity test.On the other hand, our estimates of appendix Aeffectively bound from above theprobability that $\det A \neq 0$ provided that$\det A=0 \bmod m_1$ for a random prime $m_1$ in a fixed interval.These observations  turn the computation of$r_1$ into a randomized singularity test for the matrix $A$.\end{remark}\section{\bf Computing absolutely smallerdeterminants.}\label{sec9}\setcounter{equation}{0}The next algorithm extends algorithm \ref{alg8.1}.Under the assumption that\begin{equation}\label{equ9.1}|\det A| < M-m_k,\end{equation}it tests if\begin{equation}\label{equ9.2}- -m_k \leq \det A  < m_k\end{equation}and, if so, computes $\det A$.The newalgorithm shares its stage 1 with algorithms \ref{alg7.1}and \ref{alg8.1}; in addition, it only performs subtractions of$k$ pairs of nonnegative integers bounded from above by $m_i$, for$i=1, \ldots, k$, and $2k-2$ integer distinctnesstests (each of which tests whether two given integers are distinct or not)for pairs of integers ranging from $-m_1$ to $m_1 -1$[cf.  (\ref{equ7.2})].\begin{algorithm}\label{alg9.1}, testing a range for $|\det A|$ andcomputation (modulo an integer) of absolutely smaller determinants.{\bf Input,} as for algorithm \ref{alg7.1}, assuming (\ref{equ9.1}).{\bf Output:} $\det A$ if (\ref{equ9.2}) holds,the word OVERSIZED otherwise.{\bf Computations.}\begin{enumerate}\itemas in algorithm \ref{alg7.1}.\itemIf\begin{equation}\label{equ9.3}r_i=r, \  i=1, \ldots, k,\end{equation}for some integer $r$, then output $r$.\itemOtherwise, compute $r_i-m_i, \  i=1, \ldots, k$.\itemIf\begin{equation}\label{equ9.4}r_i-m_i=r_- , \  i=1, \ldots, k\end{equation}for some integer $r_-$,then output $r_-$. Otherwise, output OVERSIZED.\end{enumerate}\end{algorithm}Let us show correctness of the algorithm.If $ \  0 \leq \det A <m_k$, then we have$$\det A = (\det A) \bmod m_i=r_i=r, \  \i=1, \ldots, k$$ [cf. (\ref{equ7.2})], so that the value $\det A$is correctly computed at stage 2.If $-m_k \leq \det A<0$, then$$\det A = ((\det A) \bmod m_i)-m_i=r_i-m_i=r_-, \  \i=1, \ldots, k, $$so that the value $\det A$ is correctly computed atstage 4.Let us next show that neither (\ref{equ9.3}) nor (\ref{equ9.4}) holdunder (\ref{equ9.1}) if\begin{equation}\label{equ9.5}\det A \geq m_k.\end{equation}Let us write\begin{equation}\label{equ9.6}\det A = r_i + s_im_i, \   \   i=1, \ldots, k, \end{equation}where  $0 \leq r_i=(\det A) \bmod m_i < m_i, \0 \leq s_i < M_i,\  i=1, \ldots, k$,\begin{equation}\label{equ9.7}  0< s_k < M_k -1,\end{equation}and all $s_i$ are integers [cf. (\ref{equ7.5}), (\ref{equ9.1}), and(\ref{equ9.5})]. Then we have\begin{equation}\label{equ9.8}r_i+s_im_i = r_k +s_km_k, \   i=1, \ldots, k-1.\end{equation}If (\ref{equ9.3}) holds, then $r_i=r_k=r, \   i=1, \ldots, k-1$,and (\ref{equ9.8}) implies that\begin{equation}\label{equ9.9}s_im_i=s_km_k, \  {\rm for\  all \   } i.\end{equation}Since $m_i$ and $m_k$ are relatively prime, (\ref{equ9.9})implies that $m_i$ divides $s_k$, for all $i<k$.Therefore, the integer $M_k=m_1\ldots m_{k-1}={\rm lcm} (m_1, \ldots, m_k)$divides $s_k$,since all $m_i$ are pairwise relatively prime. This contradicts(\ref{equ9.7}).Therefore, (\ref{equ9.3}) cannot hold under (\ref{equ9.5}).If (\ref{equ9.4}) holds, we similarly deduce from  (\ref{equ9.8}) that\begin{equation}\label{equ9.10}(s_i +1)m_i = (s_k+1)m_k, \   i=1, \ldots, k-1,\end{equation}which implies that $m_i$ divides $s_k+1$ for $i=1, \ldots, k-1$.Therefore, $M_k$ divides $s_k+1$, since all $m_i$ arepairwise relatively prime. Then again, we arrive at a contradiction to(\ref{equ9.7}).Therefore, (\ref{equ9.4}) cannot hold under (\ref{equ9.5}), which impliescorrectness of algorithm \ref{alg9.1} under (\ref{equ9.5}).To complete the correctness proof, it remains to show thatneither (\ref{equ9.3}) nor (\ref{equ9.4}) hold under (\ref{equ9.1})if\begin{equation}\label{equ9.11}\det A < -m_k.\end{equation}In this case, the integers $s_i$ of (\ref{equ9.6}),(\ref{equ9.8})--(\ref{equ9.10}) satisfy the bounds$-M_i \leq s_i < 0,  \  i=1, \ldots, k-1$,\begin{equation}\label{equ9.12} -M_k < s_k < -1.\end{equation}[cf. (\ref{equ7.5}), (\ref{equ9.1}), and (\ref{equ9.5})].Now, if (\ref{equ9.3}) holds, then we deduce from (\ref{equ9.9}) that$m_i$ divides $-s_k$ for $i=1, \ldots, k-1$, and consequently,$M_k$ divides $-s_k$.  This contradicts (\ref{equ9.12}).Therefore, (\ref{equ9.3}) cannot hold under (\ref{equ9.11}).If (\ref{equ9.4}) holds, then we deduce from (\ref{equ9.10}) that$M_k$ divides $-1-s_k$. Then again, we arrive at a contradiction to(\ref{equ9.12}).Therefore, (\ref{equ9.4}) cannot hold under (\ref{equ9.11}).This completes the correctness proof.  \hfill  $\Box$\section{\bf Computation modulo a fixed prime of thedeterminant of a matrix.}\label{sec10}\setcounter{equation}{0}Let $p$ denote a fixed prime. Then  $(\det A) \bmod p$ can beimmediately computed as soon as we compute modulo $p$the $PLU$ factorization of $A$ [see  (\ref{equ2.2})]. This can beachieved by means of Gaussian elimination with partial pivoting,which involves at most $(2/3)n^3 +O(n^2)$ arithmeticoperations (ops) modulo $p$ and at most $(n-1)n/2$ comparisonsof integers modulo $p$ with $0$.  (The latter comparisons are neededin order to ensure that only nonzero entries are used as pivots,so as to avoid divisions by $0$.)  If, at some elimination step,all the candidate entries of the pivotcolumn equal $0 \bmod p$, then we immediatelyoutput $\det A =0 \bmod p$.If we need to compute (det $A$) mod $p$ for a recursively updatedmatrix $A$, we may rely on a modified ({\it unnormalized})$QR$-factorization of $A$, such that\begin{equation}\label{equ10.1}A=(QR) \bmod p \ , \ \ \ \ \ Q^T Q = D \bmod p \ ,\end{equation}$D$ is a diagonal matrix, $R$ is a unit upper triangular matrix(compare [C]). Such a factorization can be computed(modulo a prime $p$) by means of the Givens algorithm modifiedrepectively.\begin{remark}\label{rem10.1}Under the models of parallel computing, the cited algorithms arerelatively slow: they use order of at  least $n$ parallel steps.Faster parallel algorithms (using polylogarithmic parallel timeand order of $n^3$ ops performed modulo $p$) can befound in [P85], [P87], [KP91], [KP92], [P92], [P92a],[BP], [P96] (cf. also [W] and [KS]).\end{remark}\section{\bf Decreasing the precision of the Chinese remaindercomputations.}\label{sec11}\setcounter{equation}{0}Let us recall algorithm \ref{alg7.1} of section \ref{sec7},where we assume that $m_1$ is relatively small compared to $M$but large compared to $k$. The computation of $r_i$at stage 1 and $v_i$ and $w_i$ at stage 2 can be performed modulo $m_i$with the precision of $\lceil \log_2 m_i \rceil$ bits,$i=1, \ldots,k$, which is small relatively to theprecision required in order to represent$M,\  M_1, \ldots, \ M_k$ and $D$.  Next, we are goingto examine the computation of $M_1, \ldots, M_k$at stage 2 and $D$ at stage 3, which requires  to performrelatively fewer arithmetic operations (ops) with ahigher precision [cf. (\ref{equ7.6})].We will modify the computations so as to perform themby using a little  more ops but a lower precision.Towards this goal, let us fix an {\em integer base $b$} suchthat $\log_2 b$ exceeds the computer precision. Thenthe integers lying in the range from 0 to $b-1$ will fitthe computer precision and will be called {\em "short" integers}and {\em $b$-integers}. The integers that exceed $b-1$ and theirnegations $-q$ will be called "long" and will be associatedwith the {\em $b$-polynomials}\begin{equation}\label{equ11.1}q(x)=\sum_i q_i x^i, \   0\leq q_i < b, \  {\rm for \  all} \  i,\end{equation}and with {\em $(-b)$-polynomials} $-q(x)$, respectively, such that\begin{equation}\label{equ11.2}q(b)=q, \  \  -q(b)=-q.\end{equation}The polynomials $q(x)$ and $-q(x)$ satisfying (\ref{equ11.1}),(\ref{equ11.2}) will be called the {\em $b$-associates ofthe integers } $q$ and $-q$, respectively, as well asof any pair of polynomials $Q(x)$ and $-Q(x)$ withinteger coefficients satisfying\begin{equation}\label{equ11.3}Q(b)=q.\end{equation}The transition from a polynomial $Q(x)$ with nonnegativecoefficients to its $b$-associate, which we callthe {\em $b$-reduction} of $Q(x)$, can be performed bymeans of first substituting the $b$-associate forevery coefficient of $Q(x)$, then replacing each term of $Q(x)$ byits $b$-associate, and finally, summing all the resulting $b$-polynomialstogether, by means of the known algorithms for multiprecisionsummation of nonnegative integers (in the base $b$ arithmetic)(see [Kn], pages 251--252, [O62], [O65]).Furthermore, any polynomial $Q(x)$ with integer coefficients canbe immediately represented as the difference $Q^+(x) - Q^-(x)$ oftwo polynomials $Q^+(x)$ and  $Q^-(x)$ with nonnegativeinteger coefficients.We may {\em $b$-reduce} $Q(x)$ by $b$-reducing these two polynomials$Q^+(x)$ and  $Q^-(x)$ and subtracting their $b$-associates from each otherby means of a known algorithm for multiprecision integers (inthe base $b$ arithmetic) (see [Kn], pages 252--253).We will define the class of {\em linear operations}, consistingof additions, subtractions, multiplications by ("short") $b$-integers,and divisions by ("short") $b$-integers (the latter divisionoperation will  only be used in our algorithms \ref{alg11.2}and \ref{alg11.3}, which can be both replaced by algorithm \ref{alg11.1},at the expense of using by factor of $O(log\ k)$ more ops with $b$-integers).We will need the following simple fact:\begin{fact}\label{fac11.1}For an integer $b>1$ and a pair of $b$-reducedpolynomials $U(x)$ and $V(x)$, thesign of $U(b) - V(b)$ coincides with the sign of the leadingcoefficient of the polynomial $U(x)-V(x)$.\end{fact}%%\begin{equation}\label{equ11.1}%%q=q(b)=\sum_i q_i b^i, \   0\leq q_i < b,%%\end{equation}%%$$%%(b-1)(\sum_{i=0}^h x^i) (\sum_{i=0}^k x^i) , $$%%$(b-1) \min (h+1,k+1)$.%%$(b-1)(h+1)< b^{1+ \log_b (h+1)}$%%($h+1$ decreases to $\log_b (h+1)$%%in one step), and,%%\begin{equation}\label{equ11.2s}%%C_{b,B}=O((k+\log_b h)h), \  \%% h=\lceil \log B / \log b\rceil.%%\end{equation}Let us  assume that\begin{equation}\label{equ11.4}m_1^2 < b,\end{equation}which holds, say, for$$ b= 2^{\lceil \log_2 (m_1^2 +1) \rceil }. $$Then, clearly, $m_i$, $r_i$, and $w_i$ of(\ref{equ7.2})--(\ref{equ7.6}) are $b$-integers for all $i$,and the following algorithm computes $M, M_1, \ldots, M_k$by using only linear operations of multiplication by $b$-integers:\begin{algorithm}\label{alg11.1}  $ \  \  ${\bf Input:} integers $m_1, \ldots, m_k$ satisfying (\ref{equ7.2}).{\bf Output:} the integers $M$ of (\ref{equ7.3}), $ M_1, \ldots, M_k$of (\ref{equ7.5}).{\bf Initialization.}Define the binary tree with the root 1 (at level 0), the leaves$ M_1, \ldots, M_k$ and other nodes filled with the products of someconsecutive moduli among $m_1, \ldots, m_k$ such that $m_i$ ismissing from the parent node if and only if it is missing fromat least one of its children nodes.{\bf Computations.}Recursively in $l$ compute the products of the moduli $m_i$,associated with all the nodes of level $l$ of the tree, for$l=d, d-1, \ldots, 1$.   Output the products $M_1, \ldots, M_k$,associated with level 1. Then compute and output $M=M_1m_1$.\end{algorithm}The transition from level $l$ to level $l-1$ requires $k$ multiplicationsby the moduli $m_i$, for $l=d-1, \ldots, 2$, and $k-2$ multiplicationsfor $l=d$.In the entire algorithm, we use a total of$dk-1=k\lceil \log_2 k\rceil -1$ multiplications, each by a$b$-integer $m_i$. All other operands are in the range from 0to $M_k$, that is, their $b$-associates are polynomials ofdegrees at most $\lfloor \frac{k-1}{2} \rfloor$. The entire computation byalgorithm \ref{alg11.1} costs $O(k^2\log k)$ multiplications of pairs of$b$-integers.The next algorithm uses only $k-1$ similar multiplications by $b$-integers$m_i$ for $i=2, \ldots, k$, and $k-1$ divisions of $M$ by $m_i$, for$i=1, 2, \ldots, k-1$, which have a total cost of $O(k^2)$ ring operationsand divisions with pairs of $b^2$-integers as operands.\begin{algorithm}\label{alg11.2}  $ \  \  ${\bf Input} and {\bf output} as in algorithm \ref{alg11.1}.{\bf Computations.}\begin{enumerate}\itemSet $m^{(1)}=m_i$ and recursively compute $m^{(i)}=m^{(i-1)}m_i$,$i=2, \ldots, k$. Output $M=m^{(k)}$, $M_k=m^{(k-1)}$.\itemFor $i=2, \ldots, k$, compute and output $M_i=M/m_i$,$i=1, \ldots, k-1$.\end{enumerate} \end{algorithm}At stage 2, we apply the following simple algorithm:\begin{algorithm}\label{alg11.3} $ \   \  ${\bf Input:} a $b$-polynomial $P(x)=\sum_{i=0}^h P_ix^i$ anda $b$-integer $d>0$.{\bf Output:} an integer $0\leq \rho < d$ and a $b$-polynomial$Q(x)$ of a degree at most $h$ such that $P(b)-dQ(b)=\rho$.{\bf Computations:}\begin{enumerate}\itemCompute the quotient $q_0$ and the remainder $\rho_0$ of thedivision of $P_h$ by $d$.\itemRecursively compute the quotients $q_i$ and the remainder $\rho_i$ ofthe division of $\rho_{i-1}+bP_{h-1}$ by $d$,  for $i=1, \ldots, h$.\itemOutput the integer $\rho=\rho_h$ and the $b$-polynomial$Q(x)=\sum_{i=0}^h q_ix^{h-i}$.\end{enumerate}\end{algorithm}Corectness of the algorithm and its cost bound of $O(k)$ arithmeticoperations with $b^2$-integers are immediately verified.Now, we recall that $r_i$ and $w_i$are $b$-integers under (\ref{equ11.4}), whereasthe $b$-associates of $M_i$ are $b$-polynomials of degreesless than $(k-1)/2$ for all  $i$. It follows that, based on algorithms\ref{alg11.2} and \ref{alg11.3}, we may evaluate the$b$-associates of $M_ir_iw_i$ for all $i$ and $D^* = \sum_{i=1}^k M_ir_iw_i$only involves $O(k^2)$ operations with $b$-integers.Let us consider the remaining stage of the reduction modulo $M$ of theinteger\begin{equation}\label{equ11.5}D^*=\sum_{i=1}^k M_ir_iw_i,  \0 \leq D^* < M\sum_{i=1}^k (m_i -1)=\kappa M,\end{equation}where  $\kappa= \sum_{i=1}^k (m_i-1) \leq km_1 - k -\sum_{i=1}^{k-1}i= (m_1 -\frac{k+1}{2})k$.The latter bound, (\ref{equ7.2}) and (\ref{equ11.4}) togetherimmediately imply that $\kappa$ is a $b$-integer.Now, we may restate our remaining task as follows: given a pair of$b$-polynomials $M(x)$ and $D^*(x)$, which are $b$-associates ofthe two integers $M$ of (\ref{equ7.3}) and $D^*$ of (\ref{equ11.5}),compute a $b$-associate $D(x)$ of the integer $D^*-M\bar{Q}= D^* \bmod M$such that \begin{equation}\label{equ11.6}0\leq D^* -\bar{Q}M < M\end{equation}for an interger $\bar{Q}$.(\ref{equ11.5}) and (\ref{equ11.6}) together imply that\begin{equation}\label{equ11.7}0 \leq \bar{Q} < \kappa <b,\end{equation}so that $\bar{Q}$ is a $b$-integer, which can be found by means ofthe binary searchin the interval from 0 to $\kappa-1$.  Having $\bar{Q}$ of(\ref{equ11.6}) computed, we immediately compute the polynomial$D^*(x)-M(x)\bar{Q}$ and then compute and output its $b$-associate$D(x)$.The binary search will be simplified due to thefollowing simple observation:\begin{fact}\label{fac11.2} For a pair of $b$-polynomials $U(x)$ and$V(x)$, the sign of the value $U(b)-V(b)$ of the polynomial$U(x)-V(x)$ at $x=b$ coincides with the sign of theleading coefficient of the latter polynomial.\end{fact}\begin{algorithm}\label{alg11.4},reduction modulo an integer by means of binary search.{\bf Input:} two integers $\kappa$ and $b$, $b>\kappa >1$, anda pair of $b$-polynomials $D^*(x)$and $M(x)$ such that  $D^*=D^*(b)$ and$M=M(b)$ satisfy (\ref{equ11.5}) and (\ref{equ7.3}).{\bf Output:} the integer $D=D^* \bmod M$ [cf. (\ref{equ7.6})]represented by the associated $b$-polynomial.{\bf Initialization:}write $Q^-=0$ and compute $Q^+=2^{\lceil \log_2 (\kappa -1) \rceil}$.{\bf Computations.}\begin{enumerate}\itemCompute the integer $Q=\lfloor (Q^+ + Q^-)/2 \rfloor$.\itemCompute and $b$-reduce the polynomial $M(x)Q$, so as toarrive at the $b$-polynomial $D_Q(x)$, such that$D_Q(b)=MQ$.\itema) If the leading coefficient of the polynomial $D^*(x)-D_Q(x)$is negative, then note that\begin{equation}\label{equ11.8}D^*=D^*(b)< D_Q(b)=MQ,\end{equation}write $Q^+=Q$, and repeat the computations starting withstage 1. \\b) Otherwise, check if the leading coefficient of the polynomial$D^*(x)-D_Q(x)-M(x)$ is nonnegative, and if so, then note that\begin{equation}\label{equ11.9}D^*=D^*(b) \geq D_Q(b)+M(b)=(Q+1)M,\end{equation}write $Q^-=Q$, and repeat the computationsstarting with stage 1.\\c) Otherwise, deduce from fact \ref{fac11.1} that neither(\ref{equ11.8}) or (\ref{equ11.9})hold, note that, consequently, (\ref{equ11.6}) holds,write $\bar{Q}=Q$,  and go to stage 4.\itemCompute the polynomial $\tilde{D}(x)=D^*(x)-M(x)\bar{Q}$.[Since (\ref{equ11.6}) holds, $\tilde{D}(b)=D^* \bmod M=D$.]Then compute and output its $b$-associate polynomial $D(x)$.\end{enumerate}\end{algorithm}Correctness of algorithm \ref{alg11.4} immediately followsfrom (\ref{equ11.5}) and fact \ref{fac11.1}.The algorithm requires at most$\lceil \log_2 (\kappa-1)\rceil  \leq\lceil \log_2((m_1- \frac{k+1}{2})k-1)\rceil$ steps ofthe binary search and a single pass through stage 4.Each step of the binary search involves $O(k)$ ring operations with$b$-integers, and so is stage 4.\begin{remark}\label{rem11.1}The number of ops involved in each step of binary search, in stage 4,and, consequently, in the entire algorithm can besubstantially decreased if, by means of a specialchoice of $b$ and $M>2|\det A|+1$, we arrive at asparse polynomial $M(x)$, say, at $M(x)=x^u +1$, for anappropriate integer $u$.\end{remark}Having computed the $b$-polynomial $D(x)$, we may easily computethe $b$-associate $D_{b,2}(x)$ of $2D(x)$ and then determine the signof the leading coefficient of the polynomial$D_{b,2}(x)-M(x)$. Due to fact \ref{fac11.1}, this will give us thesign of $2D(b) - M(b)=2D-M$, which will enable us to define thesign of $\det A$, based on (\ref{equ7.1}).\begin{remark}\label{rem11.2}For many input instances, we may decrease the arithmetic costof performing stages 2 and 3 of algorithm \ref{alg11.4} byskipping or not completing some of the $b$-reductions, as long asour computations  still enable us to test if theinequalities (\ref{equ11.7}) and(\ref{equ11.8}) hold. Similar comments apply tostage 4 and to the computation of $D_{b,2}(x)$.\end{remark}Based on combining algorithms \ref{alg7.1}, \ref{alg11.2}--\ref{alg11.4},we obtain the following result:\begin{theorem}\label{the11.1}Suppose that $2k+1$ given integers $b$, $m_1, \ldots, m_k$,$r_1, \ldots, r_k$ satisfy (\ref{equ7.2}) andthe following relations:$b> \sum_{i=1}^k (m_i -1)$ and  $r_i = D \bmod \  m_k$,$i=1, \ldots, k$, for some  unknown integer $D$ satisfying$$ |D|< \prod_{i=1}^k m_i -m_k. $$Then it suffices to use $O(k \log b)$ arithmetic operationsand comparisons with $b^2$-integers as operands in order tocompute $M=m_1 \ldots m_k$ and the $b$-associate $D(x)$ of $D \bmod M$.\end{theorem}\section{$p$-adic (Newton-Hensel's) lifting for thecomputation of matrix inverses and determinants.}\label{sec12}\setcounter{equation}{0}In this section, we will recall another customary approach(alternative to one of sections 7--9) to decreasing theprecision of algebraic computations with no rounding errors.The approach relies on $p$-adic (Newton-Hensel's) lifting.We will  first show how to apply this approach to computing$\det A$ by using order of $n^4$ ops, so that the algorithmis apparently inferior to one of section \ref{sec7} for larger $n$but may be competitive for small $n$.In section \ref{sec15}, we will show how to modify this approachso as to use order of $n^3$ ops. The approach will still have adisavantage of generally involving multiplications of pairsof nonconstant $b$-polynomials rather than onlymultiplications by $b$-integers, as in section \ref{sec11}.\begin{algorithm}\label{alg12.1}, $p$-adic lifting of matrix inverses.{\bf Input:} an $n \times n$ integer matrix $A$, an integer$p>1$, the matrix $S_0=A^{-1} \bmod p$, and a natural $h$.{\bf Output:} the matrix $A^{-1} \bmod p^H$, $H=2^h$.{\bf Computations.} \\Recursively compute\begin{equation}\label{equ12.1}S_j=S_{j-1}(2I-AS_{j-1}) \bmod p^J, \   J=2^j, \ j=1, \ldots, h.\end{equation}Output $S_h$.\end{algorithm}Correctness of the algorithm immediately follows fromthe next simple fact:\begin{fact}\label{fac12.1} {\rm [MC]. }$S_j=A^{-1} \bmod p^J, \  J=2^j$, for all $j$.\end{fact}{\bf Proof.} Fact \ref{fac12.1} follows from the matrix equation$I-AS_j=(I-AS_{j-1})^2 \bmod p^J$, implied by (\ref{equ12.1}).\hfill $\Box$At the $j$-th stage, the algorithm uses $O(n^3)$ opsperformed modulo $p^J$(with a precision of at most $\lceil J \log_2 p \rceil$bits), that is, a total of $O(hn^3)$ ops.To extend $p$-adic lifting of matrix inverses to$p$-adic lifting of $(\det A) \bmod p$, we next recall theequation\begin{equation}\label{equ12.2}\det A =1/ \prod_{k=1}^n (A^{-1}_k)_{k,k}\end{equation}(cf. e.g. [BP], page 327, and our section \ref{sec15}).Here and hereafter, $A_k$ denotesthe $k\times k$ leading principal submatrix of $A$, $k=1, \ldots,n$,so that $A_n=A$, and $(W)_{k,k}$ denotes the $(k,k)$-th entryof a matrix $W$. To apply (\ref{equ12.2}), we need to have theinverses modulo $p$ of $A_k$ available for all $k$ and for a fixed (prime) integer $p$. The existence of suchinverses for all $k$ is called {\it strong nonsingularity} of$A$ modulo $p$.We will assume that, for a given matrix $A$ and for our choiceof $p$, the matrix $A$ is strongly nonsingular modulo $p$.[With a random choice of a prime $p$in a fixed interval, we have  a convenient estimate forthe probability that $A$ is strongly nonsingularmodulo $p$ (see the appendix)].Then we extend algorithm \ref{alg12.1} to lifting $\det A$ as follows:\begin{algorithm}\label{alg12.2},$p$-adic lifting of matrix determinants.{\bf Input:} an integer $p>1$, an $n \times n$ matrix $A$, thematrices $S_{0,k}=A^{-1}_k \bmod p$ (so that the matrix $A$ isstrongly nonsingular modulo $p$), and a natural $h$.{\bf Output:} $\det A \bmod p^{2H}, \   H=2^h$.{\bf Computations.}\begin{enumerate}\itemApply algorithm \ref{alg12.1} to all pairs of matrices $A_k$and $S_{0,k}$ (replacing $A$ and $S_0$ in the input),so as to compute the matrices $S_{h,k}=A^{-1}_k \bmod p^H$ for$k=1, \ldots, n$.\itemCompute the value\begin{equation}\label{equ12.3}(1/\det A) \bmod p^{2H}=\prod_{k=1}^n[S_{h,k}(2I-A_kS_{h,k})]_{k,k} \bmod p^{2H}.\end{equation}\itemCompute and output$ (\det A) \bmod p^{2H}$ as the reciprocal of $(1/\det A)$modulo $p^{2H}$.\end{enumerate}\end{algorithm}Correctness of algorithm \ref{alg12.2} immediately follows from(\ref{equ12.2}) and from correctness of algorithm \ref{alg12.1}.The overall computational cost is dominated by order of $hn^4$ opsperformed at stage 1 with the precision of at most$\lceil H\log_2 p\rceil$ bits. At stage 2, we only need $O(n^3)$ ops(in order to premultiply $A_k$ by the last row of $S_{h,k}$, topostmultiply the resulting vector by the last column of $S_{h,k}$,to subtract the resulting value from $2(S_{h,k})_{k,k}$,for $k=1, \ldots, n$, and to multiply together the results forall $k$). These operations, as well as stage 3, are performed modulo $p^{2H}$, that is, with the precision of$2H\log_2p$ bits.In the next section, we will show how to perform the lattercomputations with a lower precision.To complete the computation of $(\det A) \bmod p^{2H}$, wewill complement algorithm \ref{alg12.2} by Gaussian eliminationmodulo $p$ applied in order to compute $S_{0,k}=A_k^{-1} \bmod p$for all $k$. In the special case, where the computations requireno pivoting and no row interchange, we compute and output the$L_kU_k$ factorizations modulo $p$ of $A_k$, simultaneously forall $k$, where $L_k$ and $U_k$ are the $k\times k$ leadingprincipal submatrices of the computed factors of $A$, that is,$L$ and $U$, respectively [cf. (\ref{equ2.1})].Then $A_k^{-1}=U^{-1}_kL^{-1}_k$, and Gaussianelimination outputs $A_k^{-1} \bmod p$ for all $k$, atthe cost of performing $O(n^3)$ ops modulo $p$.In the general case, Gaussian elimination withpartial pivoting either determines that $\det A =0 \bmod p$or outputs a $PLU$ factorization of $A$ modulo $p$, where$L$ and $U^T$ denote two lower triangular matrices and $P$is a permutation matrix, whose determinant ($1$ or $-1$) is available.Then we have $LU$ factorizations of the matrix $B=B_n=P^{-1}A$and of all its leading principal submatrices,$B_1, \ldots, B_{n-1}$.By applying algorithm \ref{alg12.2} to the matrix $B$ replacing$A$ in the input, we compute $(\det B) \bmod p^{2H}$and then immediately recover$(\det A) \bmod p^{2H}= \left(((\det B) \bmod p^{2H}\det P \right)\bmod p^{2H}$.To conclude, we will estimate $H$ and $p$ sufficient in orderto satisfy the bound $p^{2H}> 2|\det A|$. Due to Hadamard's bound,\begin{equation}\label{equ12.4}|\det A| \leq (a\sqrt{n})^n, \  \    a=\max_{i,j}|a_{i,j}|,\end{equation}which holds for any matrix $A=(a_{i,j})$,  it suffices to choose $p$ and$H$ satisfying\begin{equation}\label{equ12.5}2(a\sqrt{n})^n < p^{2H}, \  \   2H> n \log_p (a\sqrt{n})+\log_p 2.\end{equation}\section{\bf Decreasing the precision of the computation of$p$-adic lifting.}\label{sec13}\setcounter{equation}{0}At the $j$-th stage of algorithm \ref{alg12.1}, we generallyuse the precision of $\lceil j \log_2 p \rceil$ bits.For larger $j$, this value may exceed the computer precision$\beta$. Let us next extend the techniques of section \ref{sec11}in order to decrease the precision of computing so as tokeep it below the fixedlevel $\beta$, denoting the computer precision.We will choose the base $b=p^G$, where$G$ is a power of 2, $G\leq H/2$, so that $G$ divides $H/2$, and\begin{equation}\label{equ13.1}G \leq  \log_2 ((b-1)^2 nH/G) < \beta.\end{equation}$\beta$ denoting the computer precision.  We associatethe entries of the matrices$A_k \bmod p^J$ and $S_{j-1, k}\bmod  p^J$, for $J=2^j$ and forall $j$ and $k$, with the $b$-polynomials in $x$(see the definition of $b$-polynomials in section \ref{sec11});such polynomials take on the values of these entries for $x=b$and have degrees at most $(J/G)-1$. Let $A_{j,k}(x)$and $S_{j-1,k}(x)$ denote the resulting matrixpolynomials such that$$A_{j,k}(b) = A_k \bmod p^J, \  \   S_{j,k}(b)=S_{j,k} \bmod p^J.$$Then, for $J\geq G$, we associate the $p$-adic liftingstep (\ref{equ12.1}) with the computation of the matrix polynomial\begin{equation}\label{equ13.2}S_{j,n}(x) =S_{j-1,n}(x) (2I -A_{j,n}(x)S_{j-1,n}(x)) \bmod x^{J/G},\end{equation}where the matrix polynomials $A_{j,n}(x)$ and $S_{j-1,n}(x)$ aregiven as input and where $j=1,2, \ldots, h$.We will perform arithmetic operations modulo $x^{J/G}$ withthe matrix polynomials of (\ref{equ13.2}).Having completed such an operation, we will apply $b$-reductionto all the entries of the output matrix polynomial, followedby a new reduction modulo $x^{J/G}$. Due to the recursiveapplication of $b$-reductions and modular reductions,the entries of all the operands are polynomials with integercoefficients in the range from $1-b$ to $b-1$. Therefore, theoutput entries, even before their $b$-reductions, are thepolynomials with coefficients in the range from $-B$ to $B$ with$B \leq 2^{\beta}$, so that the $\beta$-bit precision sufficesin all these computations, due to (\ref{equ13.1}).Similar analysis applies to the computation of the matrixpolynomials associated with the matrices $S_{j,k}$ for all $j$and $k <n$.\section{\bf Decreasing the precision of computing$\det A$ when the inverses of $A_k$ are available.} \label{sec14}\setcounter{equation}{0}Having computed the matrix polynomials$A_{h+1,k}(x)$ and $S_{h,k}(x)$ associated with thematrices $A_k \bmod p^{2H}$ and $S_{h,k}$, we nextcompute the polynomial$$D^-(x)=\prod_{k=1}^n [S_{h,k}(x) (2I-A_{h+1,k}(x)S_{h,k}(x)]_{k,k}\bmod x^{2H/G},$$associated with $(1/\det A)\bmod p^{2H}$.Then again, we assume that reduction modulo $x^{2H/G}$ andthe extension of the $b$-reduction described atstages 5--7 of algorithm \ref{alg11.1} have been repeatedly appliedafter each arithmetic operation with polynomials in the processof computing $D^-(x)$. Finally, we apply a known algorithmfor computing the reciprocal of a polynomial (see e.g. [BP],page 24) in order to compute the polynomial$\hat{D}(x)=(1/D^-(x) ) \bmod x^{2H/G}$.Then, $\hat{D}(b)=(\det A) \bmod p^{2H}$ due to (\ref{equ12.3}).%%The coefficients of $D^-(x)$ and%%$\hat{D}(x)$ lie between $1-b$ and $b-1$%%(after the application of the extended%%$b$-reduction), but we may turn $D^-(x)$%%and $\hat{D}(x)$ into $b$-reduced polynomials,%%$x=b$), \ref{alg11.1}.\section{\bf Faster $p$-adic lifting ofthe determinants}\label{sec15}\setcounter{equation}{0}Some deficiency of algorithm \ref{alg12.2} is due toits excessive work for inverting all the $k \times k$ leadingprincipal submatrices $A_k$ of$A$, which requires order of $n^4$ ops.Let us show how to avoid this stage, so as to use$O(n^3)$ ops. Due to our argument of the end ofsection \ref{sec12}, we may assume that the matrix $A$is strongly nonsingular modulo $p$. Thenblock Gauss-Jourdan elimination applied to the$2 \times 2$ block matrix\begin{equation}\label{equ15.1}A=\left( \begin{array}{lr} B & C \\  E & G \end{array} \right)\end{equation}defines the following well-known block factorization of $A$:\begin{equation}\label{equ15.2}A=\left(\begin{array}{lr}I& O\\EB^{-1}& I   \end{array}\right)\left(\begin{array}{lc}  B & O\\O & S  \end{array} \right)\left( \begin{array}{lc}I & B^{-1}C \\ O & I  \end{array} \right),\end{equation}where\begin{equation}\label{equ15.3}S=G-EB^{-1}C.\end{equation}Moreover, strong nonsingularity of $A$ modulo $p$implies strong nonsingularity modulo $p$ of thematrices $B$ and $S$.Now, we arrive at the following  divide-and-conqueralgorithm.\begin{algorithm}\label{alg15.1}computing the determinant of a strongly nonsingular matrix.{\bf Input:} an $n \times n$ strongly nonsingular matrix $A$.{\bf Output:} $\det A$.{\bf Computation.}\begin{enumerate}\itemPartition the input matrix $A$ according to (\ref{equ15.1}),where $B$ is an $\lfloor n/2\rfloor \times  \lfloor n/2\rfloor$matrix. Compute $B^{-1}$ and $S$ of (\ref{equ15.3}) .\itemCompute $\det B$ and $\det S$.\itemCompute and output\begin{equation}\label{equ15.4}\det A=(\det B) \det S.\end{equation}\end{enumerate}\end{algorithm}Correctness of the algorithm immediately followsfrom (\ref{equ15.2}).Letting $C_D(k)$ and $C_I(k)$ denote the arithmeticcost of computing the determinant and the inverse ofany $k\times k$ strongly nonsingularmatrix, respectively, we obtain from (\ref{equ15.3})and (\ref{equ15.4})  that\begin{equation}\label{equ15.5}C_D(n) \leq C_D(\lfloor n/2 \rfloor) + C_D(\lceil n/2 \rceil) +C_I(\lfloor n/2 \rfloor) .\end{equation}The latter bounds also hold for solving the samecomputational problems modulo $p^{2H}$ for a fixed prime $p$and a positive integer $H$, assuming strong nonsingularity of$A$ modulo $p$.We may compute modulo $p^H$ the inverse  of a $k \times k$ matrixby using$C_I(k)=O(k^3 \log H)$ ops, by means of algorithm \ref{alg12.1}.Substitute this bound into (\ref{equ15.5}), and obtainrecursively that $C_D(n) = O(n^3 \log H)$,that is, we deduce the latter boundbased on recursive application of the factorization(\ref{equ15.2}) to the matrices $B$ and $S$ and oncomputing the inverses modulo $p^{2H}$ bymeans of algorithm \ref{alg12.1} [see the bound(\ref{equ12.5}) on $H$].Our techniques of section \ref{sec13} for decreasing theprecision of the computation of $(\det A) \bmod p^{2H}$are immediately extended to the case of algorithm \ref{alg15.1}.\section{\bf Discussion.  Decreasing the precision of thestraightforward algorithms.}\label{sec16}Further theoretical and expermental study may revealsome competitive approaches to determining the signof $\det A$. In particular,we may extend the approach of sections \ref{sec11}, \ref{sec12},and \ref{sec14} to any algorithm that outputs $\det A$ orits sign in the result of a sequence of ringoperations,  +, --, and {\large  $^ {_*}$}. For instance,we may do this for the algorithms of [E67] and [B68] (thoughtheir $b$-associates generally involve multiplications ofnonconstant $b$-polynomials) as well as for thestraightforward algorithms based on the decompositions$$ \det A=\sum_{i=1}^n (-1)^{i-1} a_{i,1} \det A_{i,1} \eqno (16.1) $$\hfill $\Box$or $$ \det A=\sum_{i=1}^n (-1)^{i-1} a_{1,j} \det A_{1,j} \eqno(16.2) $$\hfill $\Box$and on similar recursive decompositions of $\det A_{1,j}$ and/or$\det A_{i,1}$.Here, $A_{i,j}$ denotes an $(n-1) \times (n-1)$ submatrixof $A$ obtainedby deleting the $i$-th row and the $j$-th columns of $A$.Every multiplication in the latter algorithms has anoperand $a_{i,j}$ for some pair $(i,j)$. By choosing$b>\max_{i,j} |a_{i,j}|$, we may ensure thatthe $b$-associates of these operands are $b$-integers,so that the $b$-associates of the algorithms onlyinvolve linear operations.Hadamard's bound (\ref{equ12.4}), which holds for any$n \times n$  matrix $A$, enables us to bound similarythe absolute values of alloperands of these straightforward algorithms, so thatthe $b$-associates of these operands are polynomials of degreesless than $(n+1) \log_b (a\sqrt{n})$. These observations suggestthat the algorithms are competitive for small $n$.For large $n$, however, the algorithms use prohibitivelymany ops. Indeed, denoting the number of ops by $C_n$,we have $C_n = (2n +1) + nC_{n-1}$.  By extending thisrelation recursively, to $C_{n-1}, C_{n-2}, \ldots$,we obtain  the expression$1+C_n= n! \sum_{i=1}^n 1/i!$, which converges to$1+C_n=en!, \  \  e=2.7182\ldots$, as $n \longrightarrow \infty$.\vspace{.1in}\noindent{\bf Appendix A:Estimating the probability of singularity.}We will estimate the probability that det$A \neq 0$ provided that (det $A$) mod $p =0$ fora random prime $p$ and for a fixed integer matrix $A$.{\bf Fact A.1} [IR]. {\it Let $f(n) > 0$ for all natural $n$,$\displaystyle{ \lim_{n \rightarrow \infty} f(n) = \infty}$.Then there exist two positive constants $C$ and $n_0$such that the interval}$$\{ p: \ f(n)/n< p< f(n) \} \eqno (A.1)$${\it contains at least $f(n)/(C \log_2 f(n))$ distinct primesprovided that $n > n_0$.}{\bf Corollary A.1} [P85], [P87].{\it For a fixed nonsingular $n \times n$ matrix$A$ filled with integers, for $f(n) \geq n^{d+1}$, $d>0$,and for a random prime $p$ chosen in the interval (A.1),under the uniform probability distribution, we have}$$q = {\rm Probability } \{ {\rm det} \ A=0\ {\rm mod} \ p\} \leq(1+d)(C \log_2 |{\rm det} \ A|) /(f(n)d) \ . \eqno (A.2) $${\bf Proof.} Let $k$ distinct primes from the interval (A.1)divide $|{\rm det} \ A|$. Then their product $P$ alsodivides $|{\rm det} \ A|$. Therefore, this product $P$satisfies the bounds $(f(n)/n)^k < P \leq |{\rm det} \ A|$.It follows that $(\det A) \bmod p =0$ for at most$k \leq (\log_2 |{\rm det} \ A|)/\log_2 (f(n)/n)$ primes in (A.1),among at least $f(n)/(C \log_2 f(n))$. Hence,$$q \leq C {(\log_2 |{\rm det} \ A|)\log_2 f(n) \over f(n )\log_2 (f(n)/n) } \ , $$and the bound (A.2) follows for $f(n) \geq n^{1+d}$. $\Box$In particular, we may ensure that $$q={\rm Probability} \{ {\rm det} \ A=0 \ {\rm mod} \ p \}\leq \Delta \ , $$for a fixed $\Delta$, by choosing $p$ in theinterval (A.1) where, say, $$f(n) \geq \max \{ n^2 ,(2C  \log |\det A |)/ \Delta \} \ . %%\eqno (A.3)$$\vspace{.1in}\noindent{\bf Appendix B: Numerical experiments for Part I.}A set of experiments has been performed for numerical computation of the determinants of $n \times n$ matrices for $2 \le n \le 10$, based on computingthe $LU$, $PLU$, $PLUP'$ and $QR$ (Householder) factorizations of $A$.The algorithms have been implemented with C++ and built as a console applicationwith Microsoft Visual C++ 4.0 compiler and linker. All numerical operationshave been performed with double precision floating point arithmetic. The double precision representation of a number uses $64$ bits: $1$ for the sign, $11$ for the exponent, and $52$ for the mantissa. Its range is $\pm1.7\times10^{308}$ with at least $15$ decimal digits of precision. The test results have been collected on a Pentium-100MHz PC, running under Windows 95's DOS session. The system pseudo-random number generator functions {\bf srand()} and {\bf rand ()}have been used to generate input matrices.All the output values representing the determinants have been compared with the errorbounds $e_d^+$ and $e_d$ estimated according to (4.10) and (6.6), respectively. We considered three classes of input matrices $A$. In the cases where $|e_d^+| \ge |{\rm det\ }A|$ or $|e_d| \ge |{\rm det\ }A|$ for the computed values of det $A$, $e_d^+$ and/or $e_d$, ``failure" has been recorded.  Otherwise, ``correct solution" has been recorded. For each of  the three classes of $n \times n$ input matrices and for each $n$the average of the estimated relative error [min (1, $|e_d^+|$) if det $A=0$, min (1, $|{{e_d^+} \over {{\rm det\ } A}}|$) otherwise] is presented. Table \ref{TR1AF} and Table \ref{TRnAF} also include actual failure count in addition to ``estimated failure" countsince the input matrices are generated with their determinants known beforehand.It can be seen that the ``actual failures" occur substantially less frequentlythan the ``estimated failures".\noindent\smallskip{\bf Test 1: Random matrices}Determinants of 10000 random matrices have been evaluated. The entries of the input matrices are random integers in the interval $(-32768, 32768)$. The estimated relative errors are presented in Table \ref{TRE}, ``estimated failure"count is presented in Table \ref{TRF}.\begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline2	&  1.112e-011	&  1.836e-010	&  2.312e-014	&  6.877e-014 \\3	&  3.215e-005	&  9.644e-010	&  4.510e-013	&  7.702e-013 \\4	&  1.138e-003	&  3.007e-004	&  2.812e-010	&  1.827e-009 \\5	&  9.501e-003	&  2.677e-004	&  6.277e-012	&  3.729e-012 \\6	&  4.001e-002	&  8.723e-004	&  4.683e-011	&  1.049e-011 \\7	&  1.220e-001	&  3.056e-003	&  1.087e-010	&  1.256e-011 \\8	&  2.624e-001	&  5.316e-003	&  3.001e-010	&  1.438e-011 \\9	&  4.369e-001	&  8.361e-003	&  5.955e-010	&  1.733e-011 \\10	&  6.158e-001	&  1.210e-002	&  1.716e-009	&  1.904e-011 \\\hline \end{tabular}\caption{Average estimated relative error on 10,000 random matrices.}\label{TRE} \end{center} \end{table}\begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline2	&  0    	&  0    	&  0    	&  0     \\3	&  0    	&  0    	&  0    	&  0     \\4	&  8    	&  3    	&  0    	&  0     \\5	&  73   	&  2    	&  0    	&  0     \\6	&  332  	&  7    	&  0    	&  0     \\7	&  1046 	&  26   	&  0    	&  0     \\8	&  2353 	&  45   	&  0    	&  0     \\9	&  3999 	&  71   	&  0    	&  0     \\10	&  5819 	&  105  	&  0    	&  0     \\\hline \end{tabular}\caption{``Estimated failure" count on 10,000 random matrices.}\label{TRF} \end{center} \end{table}\noindent\smallskip{\bf Test 2: Matrix having $\pm1$ determinants}The input matrices are composed by using the following steps:\begin {enumerate}\item Compute a pair of low and upper triangular matrices $M$ and $N$, respectively, where the non-diagonal entries of $M$ and $N$ are random integers in the interval $(-p, p)$, for $p = \sqrt {2048 \over n}$;\item Set $M_{i,i}$ and $N_{i,i}$ to $1$, where $i = 1, \cdots, n$;\item Compute $A=MN$;\item Swap a random pair of rows in matrix A;\item Repeat step 4 (the previous step) $m$ times, where $m$ is a random integer in $[0, n)$.\end {enumerate}The average estimated relative errors are presented in Table \ref{TR1E}.The ``estimated failure" count recorded in 10,000 tests is presented in Table \ref{TR1F}.The exact values of the determinants are known to equal $(-1)^m$ in this case. If an algorithm outputs det $A$ satisfying $| {\rm det\ } A - (-1)^m| \ge 1$,an ``actual failure" is recorded for the test and the count is presented in Table \ref{TR1AF}. \begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline2		&  1.841e-010	&  1.154e-012	&  1.154e-012	&  3.138e-010 \\3		&  3.563e-004	&  1.440e-006	&  2.318e-007	&  1.129e-006 \\4		&  6.155e-002	&  3.287e-003	&  3.194e-004	&  1.528e-003 \\5		&  6.633e-001	&  3.270e-001	&  2.213e-001	&  6.315e-001 \\6 - 10	&  $>9.5e-001$	&  $>9.5e-001$	&  $>9.5e-001$	&  $>9.5e-001$ \\\hline \end{tabular}\caption{Average estimated relative error on random matrices having $\pm1$ determinants.}\label{TR1E} \end{center} \end{table}\begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline2		&  0    	&  0    	&  0    	&  0     \\3		&  1    	&  0    	&  0    	&  0     \\4		&  431  	&  16   	&  0    	&  0     \\5		&  5572 	&  1840 	&  816  	&  3986  \\6 - 10	&  $>9500$	&  $>9500$	&  $>9500$	&  $>9500$ \\\hline \end{tabular}\caption{``Estimated failure" count on 10,000 random matrices having $\pm1$ determinants.}\label{TR1F} \end{center} \end{table}\begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline$\le4$&0	&0	&0	&0 \\5	&10	&0	&0	&0 \\6	&25	&0	&0	&0 \\7	&46	&0	&0	&0 \\8	&78	&0	&0	&0 \\9	&76	&1	&5	&5 \\10	&141	&56	&114	&144 \\\hline \end{tabular}\caption{Actual failure count in 10,000 tests on random matrices having $\pm1$ determinants.}\label{TR1AF} \end{center} \end{table}\noindent\smallskip{\bf Test 3: Matrix having small determinants}The input matrices are composed by using the following steps:\begin {enumerate}\item Compute a pair of low and upper triangular matrices $M$ and $N$, respectively, where the non-zero entries of $M$ and $N$ are random integers in the interval $(-p, p)$, for $p = \sqrt {512 \over n}$;\item Compute $A=MN$ and $D = \prod_{i = 1}^nM_{i,i}N_{i,i}$;\item Swap a random pair of rows in matrix A;\item Repeat step 3 (the previous step) $m$ times, where $m$ is a random integer in $[0, n)$.\end {enumerate}The average estimated relative errors are presented in Table \ref{TRnE}.The ``estimated failure" count recorded in 10,000 tests is presented in Table \ref{TRnF}.The exact values of the determinants are known to equal $(-1)^mD$ in this case. If an algorithm outputs det $A$ satisfying $| {{\rm det\ } A \over D} - (-1)^m| \ge 1$,an ``actual failure" is recorded for the test and the count is presented in Table \ref{TRnAF}. \begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline2	&  7.207e-012	&  1.389e-012	&  1.645e-012	&  7.705e-012 \\3	&  6.776e-007	&  2.132e-007	&  2.261e-009	&  9.055e-009 \\4	&  5.100e-003	&  5.500e-004	&  9.923e-007	&  3.878e-006 \\5	&  1.005e-001	&  3.859e-003	&  3.315e-004	&  1.444e-003 \\6	&  4.015e-001	&  1.229e-001	&  5.013e-002	&  1.668e-001 \\7	&  6.864e-001	&  5.414e-001	&  4.761e-001	&  5.965e-001 \\8	&  8.367e-001	&  6.601e-001	&  6.566e-001	&  6.593e-001 \\9	&  9.381e-001	&  7.607e-001	&  7.510e-001	&  7.582e-001 \\10	&  9.844e-001	&  8.243e-001	&  8.010e-001	&  8.203e-001 \\\hline \end{tabular}\caption{Average estimated relative error on random matrices having small determinants.}\label{TRnE} \end{center} \end{table}\begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline2	&  0    	&  0    	&  0    	&  0     \\3	&  0    	&  0    	&  0    	&  0     \\4	&  43   	&  6    	&  0    	&  0     \\5	&  918  	&  24   	&  0    	&  0     \\6	&  3839 	&  937  	&  235  	&  1252  \\7	&  6709 	&  5360 	&  4640 	&  5974  \\8	&  8218 	&  6584 	&  6561 	&  6577  \\9	&  9295 	&  7553 	&  7491 	&  7530  \\10	&  9802 	&  8111 	&  7935 	&  8076  \\\hline \end{tabular}\caption{``Estimated failure" count on 10,000 random matrices having small determinants.}\label{TRnF} \end{center} \end{table}\begin{table} \begin{center} \begin{tabular}{|c|r|r|r|r|} \hlineMatrix Size & $LU$  & $PLU$ & $PLUP'$ & $QR$ \\\hline$\le5$& 0	&	0&	0&	0 \\6	& 18	&	0&	0&	0 \\7	& 49	&	0&	0&	0 \\8	& 53	&	0&	0&	0 \\9	& 55	&	0&	0&	0 \\10	& 71	&	0&	0&	0 \\\hline \end{tabular}\caption{Actual failure count in 10,000 tests on random matrices having small determinants.}\label{TRnAF} \end{center} \end{table}\newpage\vspace{.1in}\begin{center}{\bf REFERENCES}\end{center}\begin{enumerate}\item[[ABDPY]]F. Avnaim, J.-D. Boissonnat, O. Devillers, F. P. Preparata,M. Yvinec, Evaluating Signs of Determinants UsingSingle-Precision Arithmetic, to appear in {\bf Algorithmica.}\item[[AF]] D. Avis, K. Fukuda, A Pivoting Algorithm forConvex Hulls and Vertex Enumeration of Arrangements andPolyhedra, {\bf Discrete Comput. Geometry}, 8, 295--313, 1992.\item [[AHU]] A. V. Aho, J. E. Hopcroft, J. D. Ullman,{\bf The Design and Analysis of Computer Algorithms,}Addison-Wesley, Reading, Massachusetts, 1974.\item[[B67]]\AA. Bj\"{o}rck, Solving Linear Least Squares Problemsby Gram-Schmidt Orthogonalization, {\bf BIT}, 7, 1--21, 1967.\item[[B68]] E. H. Bareiss, Sylvester's Identity andMultistep Integer-Preserving Gaussian Elimination,{\bf Math. of Comp.}, 22, 565--578, 1968.\item[[BEPP]] H. Br\"onnimann, I. Z. Emiris, V. Y. Pan, S. Pion,Computing Exact Geometric Predicates Using Modular Arithmeticwith Single Precision, {\bf Proc. 13th Ann. ACM Symp. on Computational Geometry},ACM Press, New York, 1997.\item[[BM]] A. Borodin, I. Munro, {\bf The ComputationalComplexity of Algebraic and Numerical Problems}, AmericanElsevier, New York, 1975.\item[[BP]] D. Bini, V.Y. Pan, {\bf Polynomial and MatrixComputations, vol. 1, Fundamental Algorithms}, Birkhaeuser,Boston, 1994.%%\item[[BrP]]%%H. Br\"onnimann, V. Y. Pan, to appear.\item[[BS]] W. Baur, V. Strassen, On the Complexity ofPartial Derivatives, {\bf Theoretical Computer Science},22, 317-330, 1983.\item[[C]] K.L. Clarkson, Safe and Effective DeterminantEvaluation, {\bf Proc. 33rd Annual IEEE Symp. onFoundations of Computer Science}, 387-395, IEEEComputer Society Press, 1992.\item[[CdB]]S. D. Conte, C. de Boor, {\bf Elementary Numerical Analysis:An Algorithmic Approach,} McGraw Hill, New York, 1980.\item[[D]] J. Demmel, in {\bf From Topology toComputation: Proc. of the Smalefest} (M.W. Hirsch,J.E. Marsden, M. Shub editors), 305--316, Springer,New York, 1993.\item[[DB]]G. Dahlquist, \AA . Bj\"{o}rk, {\bf Numerical Methods}, Prentice-Hall, EnglewoodCliffs, New Jersey, 1967.\item[[Du]] E. Durand, {\bf Solutions Num\'{e}riquesdes \'{E}quations Alg\'{e}briques.Vol. II: Syst\`{e}mes de Plusieurs \'{E}quations.Valeurs Propres des Matrices}, Masson et Cie,Paris, 1961. %MR 24 \#B1754.\item[[DLa]] M. Deza, M. Laurent, Applications of CutPolyhedra, Report LIENS-92-18, {\bf Laboratoired'Informatique, Ecole Normale Superieure}, Paris,France, 1992, and {\bf J. of Computational andApplied Math.}, 55, 1, 191--216 and 55, 2, 217--247, 1994.\item[[E67]] J. Edmonds, Systems of Distinct Representativesand Linear Algebra, {\bf J. Res. Nat. Bur. Standards},Sect. B, 71, 4, 241--245, 1967.\item[[ECS]] I.Z. Emiris, J.F. Canny, R. Seidel, EfficientPerturbations for Handling Geometric Degeneracies,{\bf Algorithmica}, to appear.\item[[ES]] R. M. Erdahl, V. H. Smith (editors),{\bf Density Matrices and Density Functionals}, Proc.of the A. John Coleman Symp., Reidel, Dordrecht, 1987.\item[[FM]]G. Forsythe, C. Moler, {\bf Computer Solution ofLinear Algebraic Systems}, Prentice-Hall, EnglewoodCliffs, New Jersey, 1967.\item[[Fo]] L. Fox, {\bf An Introduction to NumericalLinear Algebra}, Clarendon Press, Oxford, 1964.\item[[FR]] K. Fukuda, V. Rosta, Combinatorial FaceEnumeration in Convex Polytopes, {\bf ComputationalGeometry, Theory and Applications}, 4, 191--198, 1994.\item[[GGMS]] P.E. Gill, G.H. Golub, W. Murray, M.A.Saunders, Methods for Modifying Matrix Factorizations,{\bf Math. of Computation}, 28, 505-535, 1974.\item[[GL]] G.H. Golub, C.F. Van Loan, {\bf MatrixComputations}, Johns Hopkins Univ. Press, Baltimore, MD, 1989.\item [[GPS]] K.A. Gallivan, R.J. Plemmons, A.H. Sameh, ParallelAlgorithms for Dense Linear Algebra Computations,{\bf SIAM Review}, 32, 1, 54-135, 1990.\item[[IR]] K. Iceland, M. Rosen, {\bf A ClassicalIntroduction to Modern Number Theory}, Springer,Berlin, 1982.\item [[KM]]U. Kulish, W. L. Miranker, {\bf Computer Arithmetic in Theory andPractice}, Academic Press, New York, 1980.\item[[Kn]] D.E. Knuth, {\bf The Art of Computer Programming:Seminumerical Algorithms}, Additon-Wesley, Reading, Massachusetts, 1981 and 1977.\item[[KP91]]E. Kaltofen, V. Y. Pan, Processor Efficient ParallelSolution of Linear Systems over an Abstract Field,{\bf Proc. 3rd Ann. ACM Symp. on Parallel Algorithmsand Architectures}, 180--191, ACM Press, New York, 1991.\item[[KP92]]E. Kaltofen, V. Y. Pan, Processor Efficient ParallelSolution of Linear Systems II. The Positive Characteristicand Singular Cases, {\bf Proc. 33rd Ann. IEEE Symp. onFoundation of Computer Science,} 714--723, IEEE ComputerSociety Press, 1992.%%\item[[KS]] E. Kaltofen, B.D. Saunders, On Wiedemann'sMethod for Solving Sparge Linear Systems, {\bf Proc. AAECC-5,Lecture Notes in Computer Science}, 536, 29-38, Springer,Berlin, 1991.\item[[L]] S. Linnainmaa, Taylor Expansion of the AccumulatedRounding Errors, {\bf BIT}, 16, 146-160, 1976.\item[[Ma]] R. H. Macmillan, A New Method for the NumericalEvaluation of Determinants, {\bf J. Roy.Aeronaut. Soc.}, 59, 772ff, 1955.\item[[Mu]] T. Muir, {\bf The Theory of Determinants inHistorical Order of Development}, four volumes bound astwo (I, 1693--1841; II, 1841--1860; III, 1861--1880;IV, 1881--1900), Dover, New York, 1960; {\bf Contributionsto the History of Determinants}, 1900--1920, Blackie and Son,London, 1930 and 1950.\item[[MA]] R. E. M. Moore, I. O. Angell, Voronoi Polygonsand Polyhedra, {\bf J. of Computational Physics},105, 301--305, 1993.\item[[MC]]R.T. Moenck, J.H. Carter, Approximate Algorithms to DeriveExact Solutions to Systems of Linear Equations, {\bf Proc.EUROSAM, Lecture Notes in Computer Science}, 72, 65-73,Springer, Berlin, 1979.\item[[O62]]Y. P. Ofman, On the Algorithmic Complexity of Discrete Functions,{\bf Dokl Acad. Nauk SSSR}, 145, 1, 48-51, 1962. (English transl.: {\bf Soviet Physics-Doklady}, 7, 7, 589-591, 1963).\item[[O65]]Y. P. Ofman, Universal Automaton (in Russian), {\bf Moscow Math. Society (Trudy)},14, 186-199, 1965.\item[[P84]]V. Y. Pan, {\bf How to Multiply Matrices Faster}, Lecture Notes inComputer Science, 179, Springer, 1984.\item[[P85]]V.Y. Pan, Fast and Efficient Parallel Algorithms for theExact Inversion of the Integer Matrices, {\bf Proc.5th Conference on Foundation of Software Technology andTheoretical Computer Science, Lecture Notes in ComputerScience}, 206, 504-521, Springer, Berlin, 1985.\item[[P87]]V.Y. Pan, Complexity of Parallel Matrix Computations,{\bf Theoretical Computer Science}, 54, 65-85, 1987.\item[[P88]]V.Y. Pan, Computing the Determinant and the CharacteristicPolynomial of a Matrix via Solving Linear Systems of Equations,{\bf Information Processing Letters}, 28, 2, 71-75, 1988.\item[[P92]]V.Y. Pan, Parametrization of Newton's Iteration forComputations with Structured Matrices and Applications,{\bf Computers and Math. (with Applications)}, 24, 3, 61-75, 1992.\item[[P92a]]V.Y. Pan, Complexity of Computations with Matrices andPolynomials, {\bf SIAM Review}, 34, 2, 225-262, 1992.\item[[P96]]V. Y. Pan, Parallel Computation of Polynomial GCD andSome Related Computations over Abstract Fields, {\bf TheoreticalComputer Science}, 162, 2, 173--223, 1996.\item[[R52]] J. B. Rosser, A Method of Computing ExactInverses of Matrices with Integer Coefficients,{\bf J. Res. Na. Bur. Standards}, Sect. B, 49, 349--358, 1952.\item[[S77]]G.W. Stewart, Perturbation Bounds for the QR Factorization of a Matrix,{\bf SIAM J. on Numerical Analysis}, 14, 509-18, 1977.\item[[SB]]J. Stoer, R. Bulirsch, {\bf Introduction to Numerical Analysis,}Springer-Verlag, New York, 1992.%%\item[[UP]]%%S. Ursic, C. Patarra, Exact Solution%%of Systems of Linear Equations with%%Iterative Methods, {\bf SIAM J. on Alg.%%Discr. Comp.}, 4, 111-115, 1983.\item[[W]]D.W. Wiedemann, Solving Sparse Linear Equations over FiniteFields, {\bf IEEE Trans. on Information Theory},IT-32, 1, 54-62, 1986.%%\item[[Wa]] P. Wang, A $p$-adic Algorithm%%for Univariate Partial Fractions, {\bf Proc.%%ACM Symp. on Symb. and Algebraic%%Computation}, 212-217, ACM Press, New York, 1981.\item[[Wil]]J.H. Wilkinson, {\bf The Algebraic Eigenvalue Problem,}Clarendion Press, Oxford, 1965.\item[[Y]]C. Yap, Towards Exact Geometric Computation, {\bf ComputationalGeometry, Theory and Applications,} 1996, to appear.%%%\item[[Z]]%%%R.E. Zippel, Probabilistic Algorithms for Sparse%%%Polynomials, {\bf Proc. EUROSAM 79, Lecture Notes in%%%Computer Science}, 72, 216-226, Springer, Berlin, 1979.\end{enumerate}\end{document}------- End of Forwarded Message