{'fraction_non_alphanumeric': 0.044318300957059624, 'fraction_numerical': 0.02299039944882951, 'mean_word_length': 4.712036957823595, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chapter 10: Algebraic Algorithms

Our Chapter in the upcoming Volume I: Computer Science and Software Engineering of Computing Handbook (Third edition), Allen Tucker, Teo Gonzales and Jorge L. Diaz-Herrera, editors, covers Algebraic Algorithms, both symbolic and numerical, for matrix computations and root-finding for polynomials and systems of polynomials equations. We cover part of these large subjects and include basic bibliography for further study. To meet space limitation we cite books, surveys, and comprehensive articles with pointers to further references, rather than including all the original technical papers.


Introduction
Algebraic algorithms deal with numbers, vectors, matrices, polynomials, formal power series, exponential and differential polynomials, rational functions, algebraic sets, curves and surfaces. In this vast area, manipulation with matrices and polynomials is fundamental for modern computations in Sciences and Engineering. The list of the respective computational problems includes the solution of a polynomial equation and linear and polynomial systems of equations, univariate and multivariate polynomial evaluation, interpolation, factorization and decompositions, rational interpolation, computing matrix factorization and decompositions (which in turn include various triangular and orthogonal factorizations such as LU, PLU, QR, QRP, QLP, CS, LR, Cholesky factorizations and eigenvalue and singular value decompositions), computation of the matrix characteristic and minimal polynomials, determinants, Smith and Frobenius normal forms, ranks, and (generalized) inverses, univariate and multivariate polynomial resultants, Newton's polytopes, greatest common divisors, and least common multiples as well as manipulation with truncated series and algebraic sets.
Such problems can be solved by using the error-free symbolic computations with infinite precision. Computer algebra systems such as Maple and Mathematica compute the solutions based on various nontrivial computational techniques such as modular computations, the Euclidean algorithm and continuous fraction approximation, Hensel's and Newton's lifting, Chinese Remainder algorithm, elimination and resultant methods, and Gröbner bases computation. The price to achieve perfect accuracy is the substantial memory space and computer time required to support the computations.
The alternative numerical methods rely on operations with binary or decimal numbers truncated or rounded to a fixed precision. Operating with the IEEE standard floating point numbers represented with single or double precision enables much faster computations using much smaller memory but requires theoretical and/or experimental study of the impact of rounding errors on the output. The study involves forward and backward error analysis, linear and nonlinear operators, and advanced techniques from approximation and perturbation theories. Solution of some problems involves more costly computations with extended precision. The resulting algorithms support high performance libraries and packages of subroutines such as those in Matlab, NAG SMP, LA-PACK, ScaLAPACK, ARPACK, PARPACK, MPSolve, and EigenSolve.
In this chapter we cover both approaches, whose combination frequently increases their power and enables more effective computations. We focus on the algebraic algorithms in the large, popular and highly important fields of matrix computations and root-finding for univariate polynomials and systems of multivariate polynomials. We cover part of these huge subjects and include basic bibliography for further study. To meet space limitation we cite books, surveys, and comprehensive articles with pointers to further references, rather than including all the original technical papers. Our expositions in Sections 2 and 3 follow the line of the first surveys in this area in [163,168,173,174,175].
We state the complexity bounds under the random access machine (RAM) model of computation [1,96]. In most cases we assume the arithmetic model, that is we assign a unit cost to addition, subtraction, multiplication, and division of real numbers, as well as to reading or writing them into a memory location. This model is realistic for computations with a fixed (e.g., the IEEE standard single or double) precision, which fits the size of a computer word, and then the arithmetic model turns into the word model [96]. In other cases we allow working with extended precision and assume the Boolean or bit model, assigning a unit cost to every Boolean or bitwise operation. This accounts for both arithmetic operations and the length (precision) of the operands. We denote the bounds for this complexity by O B (·). We explicitly specify whether we use the arithmetic, word, or Boolean model unless this is clear from the context. We write ops for "arithmetic operations", "log" for "log 2 " unless specified otherwise, andÕ B (·) to show that we are ignoring logarithmic factors.

Matrix Computations
Matrix computations is the most popular and a highly important area of scientific and engineering computing. Most frequently they are performed numerically, with values represented using the IEEE standard single or double precision. In the chapter of this size we must omit or just barely touch on many important subjects of this field. The reader can find further material and bibliography in the surveys [163,168] and the books [6,8,21,56,62,64,103,110,178,223,228,240] and for more specific subject areas in [6,103,223,228,238,240] on eigendecomposition and SVD, [8,56,62,103,110,223,228] on other numerical matrix factorizations, [23,130] on the over-and under-determined linear systems, their least-squares solution, and various other numerical computations with singular matrices, [106] on randomized matrix computations, [114,178] on structured matrix computations, [21,103,172,211] on parallel matrix algorithms, and [43,47,66,67,96,97,115,116,169,172,202,183,227,239] on "Error-free Rational Matrix Computations", including computations over finite fields, rings, and semirings that produce solutions to linear systems of equations, matrix inverses, ranks, determinants, characteristic and minimal polynomials, and Smith and Frobenius normal forms.

Dense, Sparse and Structured Matrices. Their Storage and Multiplication by Vectors
is an nth dimensional column vector. The straightforward algorithm computes the product Av by performing (2n − 1)m ops; this is optimal for general (dense unstructured) m × n matrices, represented with their entries, but numerous applications involve structured matrices represented with much fewer than mn scalar values. A matrix is singular if its product by some vectors vanish; they form its null space.
An m × n matrix is sparse if it is filled mostly with zeros, having only φ = o(mn) nonzero entries. An important class is the matrices associated with graphs that have families of small separators [102,134]. This includes banded matrices [b i,j ] i,j with small bandwidth 2w+1 such that b i,j = 0 unless |i−j| ≤ w. A sparse matrix can be stored economically by using appropriate data structures and can be multiplied by a vector fast, in 2φ − m ops. Sparse matrices arise in many important applications, e.g., to solving ordinary and partial differential equations (ODEs and PDEs) and graph computations.
Dense structured n × n matrices are usually defined by O(n) parameters, and one can apply FFT to multiply such matrix by a vector by using O(n log n) or O(n log 2 n) ops [178]. Such matrices are omnipresent in applications in signal and image processing, coding, ODEs, PDEs, particle simulation, and Markov chains. Most popular among them are the Toeplitz matrices T = [t i,j ] m,n i,j=1 and the Hankel matrices H = [h i,j ] m,n i,j=1 where t i,j = t i+1,j+1 and h i,j = h i+1,j−1 for all i and j in the range of their definition. Each such matrix is defined by m + n − 1 entries of its first row and first or last column. Products T v and Hv can be equivalently written as polynomial products or vector convolutions; their FFT-based computation takes O((m + n) log(m + n)) ops per product [1,21,178]. Many other fundamental computations with Toeplitz and other structured matrices can be linked to polynomial computations enabling acceleration in both areas of computing [17,18,19,20,21,24,81,85,150,151,152,168,172,178,187,208,209]. Similar properties hold for Vandermonde matrices and Cauchy matrices C = [ 1 si−tj i,j ] m,n i,j=1 where s i and t j denote m + n distinct scalars.
One can extend the structures of Hankel, Bézout, Sylvester, Frobenius (companion), Vandermonde, and Cauchy matrices to more general classes of matrices by associating linear displacement operators. (See [21,178] for the details and the bibliography.) The important classes of semiseparable, quasiseparable and other rank structured m × n matrices generalize banded matrices and their inverses; they are expressed by O(m + n) parameters and can be multiplied by vectors by performing O(m + n) ops [68,232].

Matrix Multiplication, Factorization, Randomization
The straightforward algorithm computes the m × p product AB of m × n by n × p matrices by using 2mnp − mp ops, which is 2n 3 − n 2 if m = n = p. This upper bound is not sharp. Strassen decreased it to O(n 2.81 ) ops in 1969. His result was first improved in [162] and 10 times afterward, most recently by Coppersmith and Winograd in [48], Stothers in [225], and Vasilevska Williams in [234], who use Cn ω ops for ω < 2.376, ω < 2.374 and ω < 2.3727, respectively. Due to the huge overhead constants C, however, we have that Cn ω < 2n 3 only for enormous values n. The well recognized group-theoretic techniques [44] enable a distinct description of the known matrix multiplication algorithms, but so far have only supported the same upper bounds on the complexity as the preceding works. References [225] and [234] extend the algorithms given in Reference [48], which in turn combines arithmetic progression technique with the previous advanced techniques. Each technique, however, contributes to a dramatic increase of the overhead constant that makes the resulting algorithms practically noncompetitive.
The only exception is the trilinear aggregating technique of [161] (cf. [163]), which alone supports the exponent 2.7753 [128] and together with the Any Precision Approximation (APA) techniques of [163] was an indispensable ingredient of all algorithms that have beaten Strassen's exponent 2.81 of 1969. The triple product property (TPP), which is the basis of [44], may very well have a natural link to trilinear aggregating, although the descriptions available for the two approaches are distinct. For matrices of realistic sizes the numerical algorithms in [118], relying on trilinear aggregating, use about as many ops as the algorithms of Strassen of 1969 and Winograd of 1971 but need substantially less memory space and are more stable numerically.
The exponent ω of matrix multiplication is fundamental for the theory of computing because O(n ω ) or O(n ω log n) bounds the complexity of many important matrix computations such as the computation of det A, the determinant of an n × n matrix A; its inverse A −1 (where det A = 0); its characteristic polynomial c A (x) = det( xI − A ) and minimal polynomial m A (x), for a scalar variable x; the Smith and Frobenius normal forms; the rank, rank A; a submatrix of A having the maximal rank, the solution vector x = A −1 v to a nonsingular linear system of equations A x = v, and various orthogonal and triangular factorizations of the matrix A, as well as various computations with singular matrices and seemingly unrelated combinatorial and graph computations, e.g., pattern recognition or computing all pair shortest distances in a graph [21, p. 222] or its transitive closure [1]. Consequently, all these operations use O(n ω ) ops where theoretically ω < 2.3727 [1, chap.6], [21, chap. 2]. In practice, however, the solution of all these problems takes order of n 3 ops, because of the huge overhead constant C of all known algorithms that multiply n × n matrices in Cn ω ops for ω < 2.775, the overhead of the reduction to a matrix multiplication problem, the memory space requirements, and numerical stability problems [103].
Moreover, the straightforward algorithm for matrix multiplication remains the users' choice because it is highly effective on parallel and pipeline architectures [103,211]; on many computers it supersedes even the so called "superfast" algorithms, which multiply a pair of n × n structured matrices in nearly linear arithmetic time, namely, by using O(n log n) or O(n log 2 n) ops, where both input and output matrices are represented with their short generator matrices having O(n) entries [178].
Numerous important practical problems have been reduced to matrix multiplication because it is so effective. This has also motivated the development of block matrix algorithms (called level-three BLAS, which is the acronym for Basic Linear Algebra Subprograms).
Devising asymptotically fast matrix multipliers, however, had independent technical interest. E.g., trilinear aggregating was a nontrivial decomposition of the 3-dimensional tensor associated with matrix multiplication, and [161] was the first of now numerous examples where nontrivial tensor decompositions enable dramatic acceleration of important matrix computations [124,137,160].
The two basic techniques below extend matrix multiplication. Hereafter O denotes matrices filled with zeros; I is the square identity matrices, with ones on the diagonal and zeros elsewhere.
Suppose we seek the Krylov sequence or Krylov matrix [B i v] k−1 i=0 for an n × n matrix B and an n-dimensional vector v [103,104,239]; in block Krylov computations the vector v is replaced by a matrix. The straightforward algorithm uses (2n − 1)n(k − 1) ops, that is about 2n 3 for k = n. An alternative algorithm first computes the matrix powers and then the products of n×n matrices B 2 i by n×2 i matrices, for i = 0, 1, . . . , s: The last step completes the evaluation of the Krylov sequence in 2s + 1 matrix multiplications, by using O(n ω log k) ops overall.
Special techniques for parallel computation of Krylov sequences for sparse and/or structured matrices A can be found in [170]. According to these techniques, Krylov sequence is recovered from the solution to the associated linear system (I − A) x = v, which is solved fast in the case of a special matrix A.
Another basic idea of matrix algorithms is to represent the input matrix A as a block matrix and to operate with its blocks rather than with its entries. E.g., one can compute det A and A −1 by first factorizing A as a 2 × 2 block matrix, The 2 × 2 block triangular factors are readily invertible, det A = (det A 0,0 ) det S and (BCD) −1 = D −1 C −1 B −1 , and so the cited tasks for the input A are reduced to the same tasks for the half-size matrices A 0,0 and S. It remains to factorize them recursively. The northwestern blocks (such as A 0,0 ), called leading principal submatrices, must be nonsingular throughout the recursive process, but this property holds for the highly important class of symmetric positive definite matrices A = C T C, det C = 0, and can be also achieved by means of symmetrization, pivoting, or randomization [1, chap. 6], [21, chap. 2], [178, sects. 5.5 and 5.6]). Recursive application of (1) should produce the LDU factorization A = LDU where the matrices L and U T are lower triangular and D diagonal. Having this factorization computed , we can readily solve linear systems Ax i = b i for various vectors b i , by using about 2n 2 ops for each i, rather than 2 3 n 3 + O(n 2 ) in Gaussian elimination. Factorizations (including PLU, QR, QRP, QLP, CS, LR, Cholesky factorizations and eigenvalue and singular value decompositions) are the most basic tool of matrix computations (see, e.g., [223]), recently made even more powerful with randomization (see [106,186,188,193,194,195,196,198], and the bibliography therein). It is well known that random matrices tend to be nonsingular and well conditioned (see, e.g., [218]), that is they lie far from singular matrices and therefore [103,110,223] are not sensitive to rounding errors and are suitable for numerical computations. The solution x = A −1 b of a nonsingular linear system Ax = b of n equations can be obtained with a precision p out in OB(n 3 p + n 2 p out ) Boolean time for a fixed low precision p provided the matrix A is well conditioned; that accelerates Gaussian elimination by an order of magnitude for large n + p out . Recent randomization techniques in [106,186,188,193,194,195,196,198] extend this property to much larger class of linear systems and enhance the power of various other matrix computations with singular or ill conditioned matrices, e.g., their approximation by low-rank matrices, computing a basis for the null space of a singular matrix, and approximating such bases for nearby singular matrices. Similar results have been proved for rectangular and Toeplitz matrices.
We refer the reader to [106,218,99] on impressive progress achieved in many other areas of matrix computations by means of randomization techniques.

Solution of linear systems of equations
The solution of a linear system of n equations, Ax = b is the most frequent operation in scientific and engineering computations and is highly important theoretically. Gaussian elimination solves such a system by applying (2/3)n 3 + O(n 2 ) ops.
Both Gaussian elimination and (Block) Cyclic Reduction use O(nw 2 ) ops for banded linear systems with bandwidth O(w). One can solve rank structured linear systems in O(n) ops [68,232]; generalized nested dissection uses O(n 1.5 ) flops for the inputs associated with small separator families [134,169,202].
Likewise, we can dramatically accelerate Gaussian elimination for dense structured input matrices represented with their short generators, defined by the associated displacement operators. This includes Toeplitz, Hankel, Vandermonde, and Cauchy matrices as well as matrices with similar structures. The MBA divide-and-conquer "superfast" algorithm (due to papers by Morf of 1974/1980 and Bitmead and Anderson of 1980) solves nonsingular structured linear systems of n equations in O(n log 2 n) ops by applying the recursive 2 × 2 block factorization (1) and preserving matrix structure [21,178,191,206]. In the presence of rounding errors, however, Gaussian elimination, the MBA and Cyclic Reduction algorithms easily fail unless one applies pivoting, that is interchanges the equations (and sometimes unknowns) to avoid divisions by absolutely small numbers. A by-product is the factorization A = P LU or A = P LU P ′ , for lower triangular matrices L and U T and permutation matrices P and P ′ .
Pivoting, however, takes its toll. It "usually degrades the performance" [103, page 119] by interrupting the string of arithmetic computations with the foreign operations of comparisons, is not friendly to block matrix algorithms and updating input matrices, hinders parallel processing and pipelining, and tends to destroy structure and sparseness, except for the inputs that have Cauchy-like and Vandermonde-like structure. The latter exceptional classes have been extended to the inputs with structures of Toeplitz/Hankel type by means of displacement transformation [167,178]. The users welcome this numerical stabilization, even though it slows down the MBA algorithm by a factor of n/ log 2 n, that is from "superfast" to "fast", which is still by a factor of n faster than the solution for general unstructured inputs, which takes order n 3 ops.
Can we avoid pivoting in numerical algorithms with rounding for general, sparse and structured linear systems to achieve both numerical stability and superfast performance? Yes, for the important classes where the input matri- To these input classes Gaussian elimination, Cyclic Reduction, and the MBA algorithm can be safely applied with rounding and with no pivoting. For some other classes of sparse and positive definite linear systems, pivoting has been modified into nested dissection, Markowitz heuristic rule, and other techniques that preserve sparseness during the elimination yielding faster solution without causing numerical problems [62,101,134,169,202]. Can we extend these benefits to other input matrix classes? Every nonsingular linear system A x = b is equivalent to the symmetric positive definite ones A T A x = A T b and A A T y = b where x = Ay, but great caution is recommended in such symmetrizations because the condition number κ(A) = ||A|| 2 ||A −1 || 2 ≥ 1 is squared in the transition to the matrices A T A and AA T , which means growing propagation and magnification of rounding errors.
There are two superior directions. The algorithms of [195,196,199] avoid pivoting for general and structured linear systems by applying randomization. These techniques are recent, proposed in [186, Sect. 12.2], but their effectiveness has formal and experimental support.
A popular classical alternative to Gaussian elimination is the iterative solution, e.g., by means of the Conjugate Gradient and GMRES algorithms [10,103,104,233]. They compute sufficiently long Krylov sequences (defined in the previous section) and then approximate the solution with linear combinations The cost of computing the product of the matrix A or A T A by a vector is dominant, but it is small for structured and sparse matrices A. One can even call a matrix sparse or structured if and only if it can be multiplied by a vector fast.
Fast convergence to the solution is critical. It is not generally guaranteed but proved for some important classes of input matrices. The major challenge are the extension of these classes and the design of powerful methods for special input classes, notably multilevel methods (based on the algebraic multigrid) [149,140,201] and tensor decompositions [160,124], highly effective for many linear systems arising in discretization of ODEs, PDEs, and integral equations.
Preconditioning of the input matrices at a low computational cost accelerates convergence of iterations for many important classes of sparse and structured linear systems [10,104], and more recently, based on randomized preconditioning, for quite general as well as structured linear systems [186,188,193,194,195,196,198].
One can iteratively approximate the inverse or pseudo-inverse of a matrix [103, Section 5.5.4] by means of Newton's iteration . ., and so convergence is very fast unless ||I − M X 0 || ≥ 1 or is near 1. The cost of two matrix multiplications is dominant per an iteration step; this makes the computation fast on multiprocessors as well as in the case of structured matrices M and X i . See more on Newton's iteration, including the study of its initialization, convergence, and preserving displacement matrix structure, in [178, chapters 4 and 6], [203,185,200,205,189,182].

Symbolic Matrix Computations
Rational matrix computations for a rational or integer input (such as the solution of a linear system and computing the determinant of a matrix) can be performed with no errors. To decrease the computational cost, one should control the growth of the precision of computing. Some special techniques achieve this in rational Gaussian elimination [7,97]. As a more fundamental tool one can reduce the computations modulo a sufficiently large integer m to obtain the rational or integer output values z = p/q (e.g., the solution vector for a linear system) modulo m. Then we can recover z from two integers m and z mod m by applying the continued fraction approximation algorithm, in other contexts called Euclidean algorithm [96,237]. Instead we can readily obtain z = z mod m if z mod m < r or z = −m + z mod m if z mod m < r otherwise, provided we know that the integer z lies in the range [−r, r] and if m > 2r.
Computing the determinant of an integer matrix, we can choose the modulus m based on Hadamard's bound. A nonsingular linear system Ax = v can become singular after the reduction modulo a prime p but only with a low probability for a random choice of a prime p in a fixed sufficiently large interval as well as for a reasonably large power of two and a random integer matrix [206].
One can choose m = m 1 m 2 · · · m k for pairwise relatively prime integers m 1 , m 2 , . . . , m k (we call them coprimes), then compute z modulo all these coprimes, and finally recover z by applying the Chinese Remainder algorithm [1,96]. The error-free computations modulo m i require the precision of log m i bits; the cost of computing the values z mod m i for i = 1, . . . , k dominates the cost of the subsequent recovery of the value z mod m.
Alternatively one can apply p-adic (Newton-Hensel ) lifting [96]. For solving linear systems of equations and matrix inversion they can be viewed as the symbolic counterparts to iterative refinement and Newton's iteration of the previous section, both well known in numerical linear algebra [183].
Newton's lifting begins with a prime p, a larger integer k, an integer matrix M , and its inverse Q = M −1 mod p, such that I − QM mod p = 0. Then one writes X 0 = Q, recursively computes the matrices X j = 2X j−1 − X j−1 M X j−1 mod (p 2 j ), notes that I − X j M = 0 mod (p 2 j ) for j = 1, 2, . . . , k, and finally recovers the inverse matrix M −1 from X k = M −1 mod p 2 k .
Hensel's lifting begins with the same input complemented with an integer vector b. Then one writes r (0) = b, recursively computes the vectors Newton's and Hensel's lifting are particularly powerful where the input matrices M and M −1 are sparse and/or structured, e.g., Toeplitz, Hankel, Vandermonde, Cauchy. Hensel's lifting enables the solution in nearly optimal time under both Boolean and word models [183]. We can choose p being a power of two and use computations in the binary mode. Reference [69] discusses lifting for sparse linear systems.

Computing the Sign and the Value of a Determinant
The value or just the sign of det A, the determinant of a square matrix A, are required in some fundamental geometric and algebraic/geometric computations such as the computation of convex hulls, Voronoi diagrams, algebraic curves and surfaces, multivariate and univariate resultants and Newton's polytopes. Faster numerical methods are preferred as long as the correctness of the output can be certified. In the customary arithmetic filtering approach, one applies fast numerical methods as long as they work and, in the rare cases when they fail, shifts to the slower symbolic methods. For fast numerical computation of det A one can employ factorizations A = P LU P ′ (see Section 2.2) or A = QR [45,103], precondition the matrix A [186], and then certify the output sign [207].
If A is a rational or integer matrix, then the Chinese Remainder algorithm of the previous subsection is highly effective, particularly using heuristics for working modulo m for m much smaller than Hadamard's bound on | det A| [26].
Alternatively [165,166,70], one can solve linear systems Ay(i) = b(i) for random vectors b(i) and then apply Hensel's lifting to recover det A as a least common denominator of the rational components of all y(i).
Storjohann in [224] advanced randomized Newton's lifting to yield det A more directly in the optimal asymptotic Boolean time O B (n ω+1 ) for ω < 2.3727. Wiedemann in 1986, Coppersmith in 1994, and a number of their successors compute det A by extending the Lanczos and block Lanczos classical algorithms. This is particularly effective for sparse or structured matrices A and in further extension to multivariate determinants and resultants (cf. [117,85,86,180]).

Computational Complexity Issues
Approximate solution of an nth degree polynomial equation, that is the approximation of the roots z 1 , . . . , z n for given coefficients p 0 , . . . , p n , is a classical problem that has greatly influenced the development of mathematics and computational mathematics throughout four millennia, since the Sumerian times [173,174]. The problem remains highly important for the theory and practice of the present day algebraic and algebraic/geometric computation, and new root-finding algorithms appear every year [141,142,143,144].
To approximate even a single root of a monic polynomial p(x) within error bound 2 −b we must process at least (n + 1)nb/2 bits of the input coefficients p 0 , . . . , p n−1 . Indeed perturb the x-free coefficient of the polynomial (x − 6/7) n by 2 −bn . Then the root x = 6/7 jumps by 2 −b , and similarly if we perturb the coefficients p i by 2 (i−n)b for i = 1, . . . , n−1. Thus to ensure the output precision of b bits, we need an input precision of at least (n − i)b bits for each coefficient p i , i = 0, 1, . . . , n − 1. We need at least ⌈(n + 1)nb/4⌉ bitwise operations to process these bits, each operation having at most two input bits.
It can be surprising, but we can approximate all n roots within 2 −b by using bn 2 Boolean (bit) operations up to a polylogarithmic factor for b of order n log n or higher, that is we can approximate all roots about as fast as we write down the input. We achieve this by applying the divide-and-conquer algorithms in [171,173,179] (see [123,157,219] on the related works). The algorithms first compute a sufficiently wide root-free annulus A on the complex plane, whose exterior and interior contain comparable numbers of the roots, that is the same numbers up to a fixed constant factor. Then the two factors of p(x) are numerically computed, that is F (x), having all its roots in the interior of the annulus, and G(x) = p(x)/F (x), having no roots there. Then the polynomials F (x) and G(x) are recursively factorized until factorization of p(x) into the product of linear factors is computed numerically. From this factorization, approximations to all roots of p(x) are obtained. For approximation of a single root see the competitive algorithms of [177].
It is interesting that, up to polylog factors, both lower and upper bounds on the Boolean time decrease to bn [179] if we only seek the factorization of p(x), that is, if instead of the roots z j , we compute scalars a j and b j such that The isolation of the zeros of a polynomial p(x) of (2) having integer coefficients and simple zeros is the computation of n disjoint discs, each containing exactly one root of p(x). This can be a bottleneck stage of root approximation because one can contract such discs by performing a few subdivisions and then apply numerical iterations (such as Newton's) that would very rapidly approximate the isolated zeros within a required tolerance. Reference [184] yields even faster refinement by extending the techniques of [171,173,179].
The algorithms of [171,173,179] incorporate the techniques of [157,219], but advance them and support substantially smaller upper bounds on the computational complexity. In particular these algorithms decrease by a factor of n the estimates of [219, Theorems 2.1, 19.2 and 20.1] on the Boolean complexity of polynomial factorization, root approximation and root isolation.

Root-Finding via Functional Iterations
About the same record complexity estimates for root-finding would be also supported by some functional iteration algorithms if one assumes their convergence rate defined by ample empirical evidence, although never proved formally. The users accept such an evidence instead of the proof and prefer the latter algorithms because they are easy to program and have been carefully implemented; like the algorithms of [171,173,177,179] they allow tuning the precision of computing to the precision required for every output root, which is higher for clustered and multiple roots than for single isolated roots.
For approximating a single root z, the current practical champions are modifications of Newton's iteration, z(i + 1) = z(i) − a(i)p(z(i))/p ′ (z(i)), a(i) being the step-size parameter [136], Laguerre's method [94,107], and the Jenkins-Traub algorithm [112]. One can deflate the input polynomial via its numerical division by x− z to extend these algorithms to approximating a small number of other roots. If one deflates many roots, the coefficients of the remaining factor can grow large as, e.g., in the divisor of the polynomial p(x) = x 1000 + 1 that has degree 498 and shares with p(x) all its roots having positive real parts.
For the approximation of all roots, a good option is the Weierstrass-Durand-Kerner's (hereafter WDK) algorithm, defined by the recurrence , j = 1, . . . , n, l = 0, 1, . . . . (4) It has excellent empirical global convergence. Reference [209] links it to polynomial factorization and adjusts it to approximating a single root in O(n) ops per step. A customary choice of n initial approximations z j (0) to the n roots of the polynomial p(x) (see [16] for a heuristic alternative) is given by z j (0) = r t exp(2π √ −1/n) , j = 1, . . . , n. Here t > 1 is a fixed scalar and r is an upper bound on the root radius, such that all roots z j lie in the disc {x : |x| = r} on the complex plane. This holds, e.g., for For a fixed l and for all j the computation in (4) uses O(n 2 ) ops. We can use just O(n log 2 n) ops if we apply fast multipoint polynomial evaluation algorithms based of fast FFT based polynomial division [1,21,25,178,190], but then we would face numerical stability problems. As with Newton's, Laguerre's, Jenkins-Traub's algorithms and the Inverse Power iteration in [17,208], one can employ this variant of the WDK to approximate many or all roots of p(x) without deflation. Toward this goal, one can concurrently apply the algorithm at sufficiently many distinct initial points z j (0) = r t exp(2π √ −1/N ) , j = 1, . . . , N ≥ n (on a large circle for large t) or according to [16]. The work can be distributed among processors that do not need to interact with each other until they compute the roots.
See [141,142,143,144,173] and references therein on this and other effective functional iteration algorithms. Reference [16] covers MPSolve, the most effective current root-finding subroutines, based on Ehrlich-Aberth's algorithm.

Matrix Methods for Polynomial Root-Finding
By cautiously avoiding numerical problems [103,Sec.7.4.6], one can approximate the roots of p(x) as the eigenvalues of the associated (generalized) companion matrices, that is matrices having characteristic polynomial p(x). Then one can employ numerically stable methods and the excellent software available for matrix computations, such as the QR celebrated algorithm. E.g., Matlab's subroutine roots applies it to the companion matrix of a polynomial. Fortune in [93] and in his root-finding package EigenSolve (citing earlier work of 1995 by Malek and Vaillancourt) apply it to other generalized companion matrices and update them when the approximations to the roots are improved.
The algorithms of [17,18,19,181,15,231,208,197] exploit the structure of (generalized) companion matrices, e.g., where they are diagonal plus rank-one (hereafter DPR1) matrices, to accelerate the eigenvalue computations. The papers [17,208] apply and extend the Inverse Power method [103, Section 7.6.1]; they exploit matrix structure, simplify the customary use of Rayleigh quotients for updating approximate eigenvalues, and apply special preprocessing techniques. For both companion and DPR1 inputs the resulting algorithms use linear space and linear arithmetic time per iteration step, enable dramatic parallel acceleration, and deflate the input in O(n) ops; for DPR1 matrices repeated deflation can produce all n roots with no numerical problems.
The algorithms of [18,19,15,231] employ the QR algorithm, but decrease the arithmetic time per iteration step from quadratic to linear by exploiting the rank matrix structure of companion matrices. Substantial further refinement of these techniques is required to make them competitive with MPSolve. See [244] on recent progress.
The papers [181,197] advance Cardinal's polynomial root-finders of 1996, based on repeated squaring. Each squaring is reduced to performing a small number of FFTs and thus uses order n log n ops. One can weigh potential advantage of convergence to nonlinear factors of p(x), representing multiple roots or root clusters, at the price of increasing the time per step by a factor of log n versus the Inverse Power method, advanced for root-finding in [17,208].

Extension to Approximate Polynomial GCDs
Reference [176] combines polynomial root-finders with algorithms for bipartite matching to compute approximate univariate polynomial greatest common divisor (GCD) of two polynomials, that is, the GCD of the maximum degree for two polynomials of the same or smaller degrees lying in the ǫ-neighborhood of the input polynomials for a fixed positive ǫ. Approximate GCDs are required in computer vision, algebraic geometry, computer modeling, and control. For a single example, GCD defines the intersection of two algebraic curves defined by the two input polynomials, and approximate GCD does this under input perturbations of small norms. See [14] on the bibliography on approximate GCDs, but see [167,178,195] on the structured matrix algorithms involved.

Univariate Real Root Isolation and Approximation
In some algebraic and geometric computations, the input polynomial p(x) has real coefficients, and only its real roots must be approximated. One of the fastest real root-finders in the current practice is still MPSolve, which uses almost the same running time for real roots as for all complex roots. This can be quite vexing, because very frequently the real roots make up only a small fraction of all roots [77]. Recently, however, the challenge was taken in the papers [208,197], whose numerical iterations are directed to converge to real and nearly real roots. This promises acceleration by a factor of d/r where the input polynomial has d roots, of which r roots are real or nearly real. In the rest of this section we cover an alternative direction, that is real root-finding by means of isolation of the real roots of a polynomial.
We write p(x) = a d x d + · · · + a 1 x + a 0 , assume integral coefficients with the maximum bit size τ = 1 + max i≤d {lg |a i |}, and seek isolation of real roots, that is seek real line intervals with rational endpoints, each containing exactly one real root. We may seek also the root's multiplicity. We assume rational algorithms, that is, error-free algorithms that operate with rational numbers.
If all roots of p(x) are simple, then the minimal distance between them, the separation bound, is at most b = d −(d+2)/2 (d + 1) (1−d)/2 2 τ (1−d) , or roughly 2 − O(dτ ) (e.g., [147]), and we isolate real roots as soon as we approximate them within less than b/2. Effective solution algorithms rely on Continued Fractions (see below), having highly competitive implementation in synaps [153,109] and its descendant realroot, a package of mathemagix, on the Descartes' rule of signs, and the Sturm or Sturm-Habicht sequences. The bounds exceed those of Theorem 3.1, but [204] has changed this, by closing the gap. Moreover rational solvers are heavily in use, have long and respected history, and are of independent technical interest. Most popular are the subdivision algorithms, such as sturm, descartes and bernstein. By mimicking binary search, they repeatedly subdivide an initial interval that contains all real roots until every tested interval contains at most one real root. They differ in the way of counting the real roots in an interval.
The algorithm sturm (due to the work by Sturm of 1835, see [96]) is the closest to binary search; it produces isolating intervals and root multiplicities at the cost O B (d 4 τ 2 ) [63,83]; see [77] on the decrease of the expected cost to O B (rd 2 τ ).
The complexity of both algorithms descartes and bernstein is O B (d 4 τ 2 ) [72,83]. Both rely on Descartes' rule of sign, but the bernstein algorithm also employs the Bernstein basis polynomial representation. See [235,2] on the theory and history of descartes, [46,216,71,146,217] on its modern versions, and [83,156] and the references therein on the bernstein algorithm.
The Continued Fraction algorithm, cf, computes the continued fraction expansions of the real roots of the polynomial. The first formulation of the algorithm is due to Vincent. By Vincent's theorem repeated transforms x → c + 1 x eventually yield a polynomial with zero or one sign variation and thus (by Descartes' rule) with zero or resp. one real root in (0, ∞). In the latter case the inverse transformation computes an isolating interval. Moreover, the c's in the transform correspond to the partial quotients of the continued fraction expansion of the real root. Variants differ in the way they compute the partial quotients.
Recent algorithms control the growth of coefficient bit-size and decrease the bit-complexity from exponential (of Vincent) to O B (d 3 τ ) expected and O B (d 4 τ 2 ) worst-case bit complexity. See [145,220,229,230] and the references therein on these results, history and variants of CP algorithms.

Systems of Nonlinear Equations
Given a system {p 1 (x 1 , . . . , x n ), . . . , p r (x 1 , . . . , x n )} of nonlinear polynomials with rational coefficients, the n-tuple of complex numbers (a 1 , . . . , a n ) is a solution of the system if p i (a 1 , . . . , a n ) = 0, 1 ≤ i ≤ r. Each p i (x 1 , . . . , x n ) is said to be an element of Q[x 1 , . . . , x n ], the ring of polynomials in x 1 , . . . , x n over the field of rational numbers. In this section, we explore the problem of solving a well-constrained system of nonlinear equations, namely when r = n, which is the typical case in applications. We also indicate how an initial phase of exact algebraic computation leads to certain numerical methods that can approximate all solutions; the interaction of symbolic and numeric computation is currently an active domain of research, e.g. [22,82,125]. We provide an overview and cite references to different symbolic techniques used for solving systems of algebraic (polynomial) equations. In particular, we describe methods involving resultant and Gröbner basis computations.
Resultants, as explained below, formally express the solvability of algebraic systems with r = n + 1; solving a well-constrained system reduces to a resultant computation as illustrated in the sequel. The Sylvester resultant method is the technique most frequently utilized for determining a common root of two polynomial equations in one variable. However, using the Sylvester method successively to solve a system of multivariate polynomials proves to be inefficient.
It is more efficient to eliminate n variables together from n + 1 polynomials, thus, leading to the notion of the multivariate resultant. The three most commonly used multivariate resultant matrix formulations are those named after Sylvester or Macaulay [36,38,135], those named after Bézout or Dixon [33,60,121], or the hybrid formulation [57,113,122]. Extending the Sylvester-Macaulay type, we shall emphasize also sparse resultant formulations [37,98,226]. For a unified treatment, see [81].
The theory of Gröbner bases provides powerful tools for performing computations in multivariate polynomial rings. Formulating the problem of solving systems of polynomial equations in terms of polynomial ideals, we will see that a Gröbner basis can be computed from the input polynomial set, thus, allowing for a form of back substitution in order to compute the common roots.
Although not discussed, it should be noted that the characteristic set algorithm can be utilized for solving polynomial systems. Although introduced for studying algebraic differential equations [214], the method was converted to ordinary polynomial rings when developing an effective method for automatic theorem proving [242]. Given a polynomial system P , the characteristic set algorithm computes a new system in triangular form, such that the set of common roots of P is equivalent to the set of roots of the triangular system [120]. Triangular systems have k 1 polynomials in a specific variable, k 2 polynomials in this and one more variable, k 3 polynomials in these two and one more variable, and so on, for a total number of k 1 + · · · + k n polynomials.

Resultant of Univariate Systems
The question of whether two polynomials f ( have a common root leads to a condition that has to be satisfied by the coefficients of f, g. Using a derivation of this condition due to Euler, the Sylvester matrix of f, g (which is of dimension m+n) can be formulated. The vanishing of the determinant of the Sylvester matrix, known as the Sylvester resultant, is a necessary and sufficient condition for f, g to have common roots in the algebraic closure of the coefficient ring.
As a running example let us consider the following bivariate system [131]: Without loss of generality, the roots of the Sylvester resultant of f and g treated as polynomials in y, whose coefficients are polynomials in x, are the x-coordinates of the common roots of f, g. More specifically, the Sylvester resultant with respect to y is given by the following determinant: An alternative matrix of order max{m, n}, named after Bézout, yields the same determinant.
The roots of the Sylvester determinant are {−3, 0, 1}. For each x value, one can substitute the x value back into the original polynomials yielding the solutions (−3, 1), (0, 1), (1, −1). More practically, one can use the Sylvester matrix to reduce system solving to the computation of eigenvalues and eigenvectors as explained in "Polynomial System Solving by Using Resultants".
The Sylvester formulations has led to a subresultant theory, which produced an efficient algorithm for computing the GCD of univariate polynomials and their resultant, while controlling intermediate expression swell [213,133]. Subresultant theory has been generalized to several variables, e.g. [32,53].

Resultants of Multivariate Systems
The solvability of a set of nonlinear multivariate polynomials is determined by the vanishing of a generalization of the resultant of two univariate polynomials. We examine two generalizations: the classical and the sparse resultants. Both generalize the determinant of n + 1 linear polynomials in n variables.
The classical resultant of a system of n + 1 polynomials with symbolic coefficients in n variables vanishes exactly when there exists a common solution in the projective space over the algebraic closure of the coefficient ring [50]. The sparse (or toric) resultant characterizes solvability of the same overconstrained system over a smaller space, which coincides with affine space under certain genericity conditions [51,98,226]. The main algorithmic question is to construct a matrix whose determinant is the resultant or a nontrivial multiple of it.
Cayley, and later Dixon, generalized Bézout's method to a set {p 1 (x 1 , . . . , x n ) , . . . , p n+1 (x 1 , . . . , x n )} of n + 1 polynomials in n variables. The vanishing of the determinant of the Bézout-Dixon matrix is a necessary and sufficient condition for the polynomials to have a nontrivial projective common root, and also a necessary condition for the existence of an affine common root [33,60,81,121]. A nontrivial resultant multiple, known as the projection operator, can be extracted via a method discussed in [41, thm. 3.3.4]. This article, along with [73], explain the correlation between residue theory and the Bézout-Dixon matrix; the former leads to an alternative approach for studying and approximating all common solutions. Macaulay [135] constructed a matrix whose determinant is a multiple of the classical resultant; he stated his approach for a well-constrained system of n homogeneous polynomials in n variables. The Macaulay matrix simultaneously generalizes the Sylvester matrix and the coefficient matrix of a system of linear equations. Like the Dixon formulation, the Macaulay determinant is a multiple of the resultant. Macaulay, however, proved that a certain minor of his matrix divides the matrix determinant to yield the exact resultant in the case of generic coefficients. To address arbitrary coefficients, Canny [36] proposed a general method that perturbs any polynomial system and extracts a nontrivial projection operator from Macaulay's construction.
By exploiting the structure of polynomial systems by means of sparse elimination theory, a matrix formula for computing the sparse resultant of n+1 poly-nomials in n variables was given in [37] and consequently improved in [40,76]. Like the Macaulay and Dixon matrices, the determinant of the sparse resultant matrix, also known as Newton matrix, only yields a projection operation. However, in certain cases of bivariate and multihomogeneous systems, determinantal formulae for the sparse resultant have been derived [57,80,122]. To address degeneracy issues, Canny's perturbation has been extended in the sparse context [54]. D'Andrea [52] extended Macaulay's rational formula for the resultant to the sparse setting, thus defining the sparse resultant as the quotient of two determinants; see [79] for a simplified algorithm in certain cases.
Here, sparsity means that only certain monomials in each of the n + 1 polynomials have nonzero coefficients. Sparsity is measured in geometric terms, namely, by the Newton polytope of the polynomial, which is the convex hull of the exponent vectors corresponding to nonzero coefficients. The mixed volume of the Newton polytopes of n polynomials in n variables is defined as an integer-valued function that bounds the number of toric common roots of these polynomials [13]. This remarkable bound is the cornerstone of sparse elimination theory. The mixed volume bound is significantly smaller than the classical Bézout bound for polynomials with small Newton polytopes but they coincide for polynomials whose Newton polytope is the unit simplex multiplied by the polynomial's total degree. Since these bounds also determine the degree of the sparse and classical resultants, respectively, the latter has larger degree for sparse polynomials. Last, but not least, the classical resultant can identically vanish over sparse systems, whereas the sparse resultant can still yield the desired information about their common roots [51].

Polynomial System Solving by Using Resultants
Suppose we are asked to find the common roots of a set of n polynomials in n variables {p 1 (x 1 , . . . , x n ), . . ., p n (x 1 , . . . , x n )}. By augmenting this set by a generic linear polynomial [36,51], we construct the u-resultant of a given system of polynomials. The u-resultant is named after the indeterminates u, traditionally used to represent the generic coefficients of the additional linear polynomial. The u-resultant factors into linear factors over the complex numbers, providing the common roots of the given polynomials equations. The method relies on the properties of the multivariate resultant, and hence, can be constructed using either Macaulay's, Dixon's, or sparse formulations. An alternative approach is to hide a variable in the coefficient field [74,81,138].
Consider the previous example augmented by a generic linear form: As described in [38], the following (transposed) Macaulay matrix M corre-sponds to the u-resultant of the above system of polynomials: It should be noted that corresponds to the affine solutions (1, −1), (−3, 1), (0, 1), whereas one solution at infinity corresponds to the last factor. Resultant matrices can also reduce polynomial system solving to a regular or generalized eigenproblem (cf. "Matrix Eigenvalues and Singular Values Problems"), thus, transforming the nonlinear question to a problem in linear algebra. This is a classical technique that enables us to numerically approximate all solutions [4,39,41,74,81]. For demonstration, consider the previous system and its resultant matrix M . The matrix rows are indexed by the following row vector of monomials in the eliminated variables: v = x 3 , x 2 y, x 2 , xy 2 , xy, x, y 3 , y 2 , y, 1 .
Vector vM expresses the polynomials indexing the columns of M , which are multiples of the three input polynomials by various monomials. Let us specialize variables u and v to random values. Then M contains a single variable w and is denoted M (w). Solving the linear system vM (w) = 0 in vector v and in scalar w is a generalized eigenproblem, since M (w) can be represented as M 0 + wM 1 , where M 0 and M 1 have numeric entries. If, moreover, M 1 is invertible, we arrive at the following eigenproblem: For every solution (a, b) of the original system, there is a vector v among the computed eigenvectors, which we evaluate at x = a, y = b and from which the solution can be recovered by division [74]. As for the eigenvalues, they correspond to the values of w at the solutions; see [75] on numerical issues, and an implementation. An alternative method for approximating or isolating all real roots of the system is to use the so-called Rational Univariate Representation (RUR) of algebraic numbers [35,215]. This allows us to express each root coordinate as the value of a univariate polynomial, evaluated over an algebraic number, which is specified as a solution of a single polynomial equation. All polynomials involved in this approach are derived from the resultant.
The resultant matrices are sparse and have quasi Toeplitz/Hankel structure (also called multilevel Toeplitz/Hankel structure), which enables their fast multiplication by vectors. By combining the latter property with various advanced nontrivial methods of multivariate polynomial root-finding, substantial acceleration of the construction and computation of the resultant matrices and approximation of the system's solutions was achieved in [24,85,86,150,151,152].
A comparison of the resultant formulations can be found, e.g., in [81,120,138]. The multivariate resultant formulations have been used for diverse applications such as algebraic and geometric reasoning [41,59,138], including separation bounds for the isolated roots of arbitrary polynomial systems [84], robot kinematics [55,212,138], and nonlinear computational geometry, computeraided geometric design and, in particular, implicitization [32,42,78,87,111].
A term t = x e1 1 x e2 2 . . . x en n of a polynomial is a product of powers with deg(t) = e 1 + · · · + e n . In order to add needed structure to the polynomial ring we will require that the terms in a polynomial be ordered in an admissible fashion [50,97]. Two of the most common admissible orderings are the lexicographic order (≺ l ), where terms are ordered as in a dictionary, and the degree order (≺ d ), where terms are first compared by their degrees with equal degree terms compared lexicographically. A variation to the lexicographic order is the reverse lexicographic order, where the lexicographic order is reversed.
Much like a polynomial remainder process, the process of polynomial reduction involves subtracting a multiple of one polynomial from another to obtain a smaller degree result [50,105,127]. A polynomial g is said to be reducible with respect to a set P = {p 1 , . . . , p r } of polynomials if it can be reduced by one or more polynomials in P . When g is no longer reducible by the polynomials in P , we say that g is reduced or is a normal form with respect to P .
For an arbitrary set of basis polynomials, it is possible that different reduction sequences applied to a given polynomial g could reduce to different normal forms. A basis G ⊆ Q[x 1 , . . . , x n ] is a Gröbner basis if and only if every polynomial in Q[x 1 , . . . , x n ] has a unique normal form with respect to G. Buchberger [27,28,29] showed that every basis for an ideal (p 1 , . . . , p r ) in Q[x 1 , . . . , x n ] can be converted into a Gröbner basis {p * 1 , . . . , p * s } = GB(p 1 , . . . , p r ), concomitantly designing an algorithm that transforms an arbitrary ideal basis into a Gröbner basis. Another characteristic of Gröbner bases is that by using the above mentioned reduction process we have g ∈ (p 1 . . . , p r ) ⇐⇒ g mod (p * 1 , . . . , p * s ) = 0 .
Further, by using the Nullstellensatz it can be shown that p 1 . . . , p r viewed as a system of algebraic equations is solvable if and only if 1 ∈ GB(p 1 , . . . , p r ). Depending on which admissible term ordering is used in the Gröbner bases construction, an ideal can have different Gröbner bases. However, an ideal cannot have different (reduced) Gröbner bases for the same term ordering. Any system of polynomial equations can be solved using a lexicographic Gröbner basis for the ideal generated by the given polynomials. It has been observed, however, that Gröbner bases, more specifically lexicographic Gröbner bases, are hard to compute [139]. In the case of zero-dimensional ideals, those whose varieties have only isolated points, a change of basis algorithm was outlined in [90], which can be utilized for solving: one computes a Gröbner basis for the ideal generated by a system of polynomials under a degree ordering. The so-called change of basis algorithm can then be applied to the degree ordered Gröbner basis to obtain a Gröbner basis under a lexicographic ordering. Significant progress has been achieved in the algorithmic realm by Faugère [88,89].
Another way to finding all common real roots is by means of RUR; see the previous section. All polynomials involved in this approach can be derived from the Gröbner basis. A rather recent development concerns the generalization of Gröbner bases to border bases, which contain all information required for system solving but can be computed faster and seem to be numerically more stable [127,154,222,155].
Turning to Lazard's example in form of a polynomial basis, one obtains (under lexicographical ordering with x≺ l y) a Gröbner basis in which the variables are triangulated such that the finitely many solutions can be computed via back substitution: The final univariate polynomial has minimal degree, whereas the polynomials used in the back substitution have total degree no larger than the number of roots. As an example, x 2 y 2 is reduced with respect to the previously computed Gröbner basis {p * 1 , p * 2 , p * 3 } = GB(p 1 , p 2 ) along two distinct reduction paths, both yielding −3x − 2y + 2 as the normal form.
There is a strong connection between lexicographic Gröbner bases and the previously mentioned resultant techniques. For some types of input polynomials, the computation of a reduced system via resultants might be much faster than the computation of a lexicographic Gröbner basis.

Research Issues and Summary
Algebraic algorithms deal with numbers, vectors, matrices, polynomials, formal power series, exponential and differential polynomials, rational functions, algebraic sets, curves and surfaces. In this vast area, manipulations with matrices and polynomials, in particular the solution of a polynomial equation and linear and polynomial systems of equations, are most fundamental in modern computations in Sciences, Engineering, and Signal and Image Processing. We reviewed the state of the art for the solution of these three tasks and gave pointers to the extensive bibliography.
Among numerous interesting and important research directions of the topics in Sections 2 and 3, we wish to cite computations with structured matrices, including their applications to polynomial root-finding, currently of growing interest, and new techniques for randomized preprocessing for matrix computations, evaluation of resultants and polynomial root-finding. Section 4 of this chapter has briefly reviewed polynomial system solving based on resultant matrices as well as Gröbner bases. Both approaches are currently active. This includes practical applications to small and medium-size systems. Efficient implementations that handle the nongeneric cases, including multiple roots and nonisolated solutions, is probably the most crucial issue today in relation to resultants. The latter are also studied in relation to a more general object, namely the discriminant of a well-constrained system, which characterizes the existence of multiple roots. Another interesting current direction is algorithmic improvement by exploiting the structure of the polynomial systems, including sparsity, or the structure of the encountered matrices, for both resultants and Gröbner bases.

Defining Terms
Characteristic polynomial: Shift an input matrix A by subtracting the identity matrix xI scaled by variable x. The determinant of the resulting matrix is the characteristic polynomial of the matrix A. Its roots coincide with the eigenvalues of the shifted matrix A − xI.
Condition number of a matrix is a scalar κ which grows large as the matrix approaches a singular matrix; then numeric inversion becomes an ill-conditioned problem. κ OUTPUT ERROR NORM ≈ INPUT ERROR NORM.
Degree order: An order on the terms in a multivariate polynomial; for two variables x and y with x ≺ y the ascending chain of terms is 1 ≺ x ≺ y ≺ x 2 ≺ xy ≺ y 2 · · ·.
Determinant: A polynomial in the entries of a square matrix whose value is invariant in adding to a row (resp. column) any linear combination of other rows (resp. columns). det(AB) = det A·det B for a pair of square matrices A and B, det B = − det A if the matrix B is obtained by interchanging a pair of adjacent rows or columns of a matrix A, det A = 0 if and only if a matrix A is invertible. Determinant of a block diagonal or block triangular matrix is the product of the diagonal blocks, and so det A = (det A 0.0 ) det S under (1). One can compute a determinant by using these properties and matrix factorizations, e.g., recursive factorization (1).
Gröbner basis: Given a term ordering, the Gröbner basis of a polynomial ideal is a generating set of this ideal, such that the (multivariate) division of any polynomial by the basis has a unique remainder.
Lexicographic order: An order on the terms in a multivariate polynomial; for two variables x and y with x ≺ y the ascending chain of terms is 1 ≺ x ≺ x 2 ≺ · · · ≺ y ≺ xy ≺ x 2 y · · · ≺ y 2 ≺ xy 2 · · ·.
Matrix eigenvector: A column vector v such that Av = λv, for a square matrix A and the associated eigenvalue λ. A generalized eigenvector v satisfies the equation Av = λBv for two square matrices A and B and the associated eigenvalue λ. Both definitions extend to row vectors that premultiply the associated matrices.
Mixed volume: An integer-valued function of n convex polytopes in n-dimensional Euclidean space. Under proper scaling, this function bounds the number of toric complex roots of a well-constrained polynomial system, where the convex polytopes are defined to be the Newton polytopes of the given polynomials.
Newton polytope: The convex hull of the exponent vectors corresponding to terms with nonzero coefficients in a given multivariate polynomial.
Resultant: A polynomial in the coefficients of a system of n polynomials with n + 1 variables, whose vanishing is the minimal necessary and sufficient condition for the existence of a solution of the system.

Separation bound:
The minimum distance between two (complex) roots of a univariate polynomial.
Singularity: A square matrix is singular if its product with some nonzero matrix is the zero matrix. Singular matrices do not have inverses.
Sparse matrix: A matrix whose zero entries are much more numerous than its nonzero entries.
Structured matrix: A matrix whose every entry can be derived by a formula depending on a smaller number of parameters, typically on O(m + n) parameters for an m×n matrix, as opposed to its mn entries. For instance, an m × n Cauchy matrix has 1 si−tj as the entry in row i and column j and is defined by m + n parameters s i and t j , i = 1, . . . , m; j = 1, . . . , n. Typically a structured matrix can be multiplied by a vector in nearly linear arithmetic time.
the graph and combinatorial computations related to linear algebra are cited in "Some Computations Related to Matrix Multiplication" and [169].
There is a lot of generic software packages for exact computation, synaps [153], a C++ open source library devoted to symbolic and numeric computations with polynomials, algebraic numbers and polynomial systems, which has been evolving into the realroot package of the open source computer algebra system mathemagix; ntl a high-performance C++ library providing data structures and algorithms for vectors, matrices, and polynomials over the integers and finite fields, and exacus [11], a C++ library for curves and surfaces that provides exact methods for solving polynomial equations. A highly efficient tool is FGb for Gröbner basis, and RS for the rational univariate representation, and real solutions of systems of polynomial equations and inequalities. Finally, Lin-Box [65] is a C++ library that provides exact high-performance implementations of linear algebra algorithms.
This chapter does not cover the area of polynomial factorization. We refer the interested reader to [96,132,158], and the bibliography therein.
The SIAM Journal on Matrix Analysis and Applications and Linear Algebra and Its Applications are specialized on Matrix Computations, Mathematics of Computation and Numerische Mathematik are leading among numerous other good journals on numerical computing.
The Journal of Symbolic Computation and the Foundations of Computational Mathematics specialize on topics in Computer Algebra, which are also covered in the Journal of Computational Complexity, the Journal of Pure and Applied Algebra and, less regularly, in the Journal of Complexity. Mathematics for Computer Science and Applicable Algebra in Engineering, Communication and Computing are currently dedicated to the subject of the chapter as well. Theoretical Computer Science has become more open to algebraic-numerical and algebraic-geometric subjects [22,34,82,125].
The annual International Symposium on Symbolic and Algebraic Computation (ISSAC) is the main conference in computer algebra; these topics are also presented at the bi-annual Conference MEGA and the newly founded SIAM conference on Applications of Algebraic Geometry. They also appear, in the annual ACM Conference on Computational Geometry, as well as at various Computer Science conferences, including SODA, FOCS, and STOC. Among many conferences on numerical computing, most comprehensive ones are organized under the auspices of SIAM and ICIAM. The International Workshop on Symbolic-Numeric Algorithms can be traced back to 1997 (SNAP in INRIA, Sophia Antipolis) and a special session in IMACS/ACA98 Conference in Prague, Czech Republic, in 1998 [175]. It restarted in Xi'an, China,