The FastTwoSum algorithm is a classical way to evaluate the rounding error that occurs when adding two numbers in finite precision arithmetic. Starting with Dekker in the early 1970s, numerous floating-point analyses have been made of this algorithm, that are aimed at identifying sufficient conditions for the error to be computed exactly and, otherwise, at quantifying the quality of the error estimate thus produced. In this paper we revisit these two aspects of FastTwoSum. We first provide new, less restrictive conditions for exactness, and show that FastTwoSum performs an error-free transform in more general situations than those found so far in the literature. Second, when exactness cannot be guaranteed we give several error analyses of the output of FastTwoSum and show that the bounds obtained are tight. In particular, this provides further insight into how the algorithm behaves when roundings other than ‘to nearest’ are used, or when the operands are reversed.

We show that the availability of fused arithmetic operators that evaluate expressions of the form ab + cd (FD2 instruction) or ab + cd + e (FD2A instruction) in floating-point arithmetic with one final rounding only would significantly facilitate many calculations that are hard to perform with high accuracy at small cost using only the traditional operations +, −, ×, ÷, √, and fused multiply-add (FMA).

We design algorithms for the correct rounding of the power function xyin the binary64 IEEE 754 format, for all rounding modes, modulo the knowledge of hardest-to-round cases. Our implementation of these algorithms largely outperforms previous correctly-rounded implementations and is not far from the efficiency of current mathematical libraries, which are not correctly-rounded. Still, we expect our algorithms can be further improved for speed. The proofs of correctness are fully detailed in the extended version [9] of this paper, with the goal to enable a formal proof of these algorithms. We hope this work will motivate the next IEEE 754 revision committee to require correct rounding for mathematical functions.

We analyze two fast and accurate algorithms recently presented by Borges for computing$x^{-1/2}$in binary floating-point arithmetic (assuming that efficient and correctly-rounded FMA and square root are available). The first algorithm is based on the Newton-Raphson iteration, and the second one uses an order-3 iteration. We give attainable relative-error bounds for these two algorithms, build counterexamples showing that in very rare cases they do not provide a correctly-rounded result, and characterize precisely when such failures happen in IEEE 754 binary32 and binary64 arithmetics. We then give a generic (i.e., precision-independent) algorithm that always returns a correctly-rounded result, and show how it can be simplified and made more efficient in the important cases of binary32 and binary64.

We review several ways to split a floating-point number, that is, to decompose it into the exact sum of two floating-point numbers of smaller precision. All the methods considered here involve only a few IEEE floating-point operations, with rounding to nearest and including possibly the fused multiply -add (FMA). Applications range from the implementation of integer functions such as round and floor to the computation of suitable scaling factors aimed, for example, at avoiding spurious underflows and overflows when implementing functions such as the hypotenuse.

We study the accuracy of classical algorithms for evaluating expressions of the form √ (a2+ b2) and c/√ (a2+ b2)in radix-2, precision-p floating-point arithmetic, assuming that the elementary arithmetic operations ±, x, /, '/ are rounded to nearest, and assuming an unbounded exponent range. Classical analyses show that the relative error is bounded by 2u+O(u2) for √ (a2+ b2), and by 3u+O(u2) for c/√ (a2+ b2), where u = 2-pis the unit roundoff. Recently, it was observed that for √ (a2+ b2) the O(u2) term is in fact not needed [1]. We show here that it is not needed either for c√ (a2+ b2). Furthermore, we show that these error bounds are asymptotically optimal. Finally, we show that both the bounds and their asymptotic optimality remain valid when an FMA instruction is used to evaluate a2+ b2.

We compute minimal bases of solutions for a general interpolation problem, which encompasses Hermite-Pade approximation and constrained multivariate interpolation, and has applications in coding theory and security. This problem asks to find univariate polynomial relations between m vectors of size σ; these relations should have small degree with respect to an input degree shift. For an arbitrary shift, we propose an algorithm for the computation of an interpolation basis in shifted Popov normal form with a cost of O~(mω-1 σ) field operations, where ω is the exponent of matrix multiplication and the notation O~(·) indicates that logarithmic terms are omitted.

Assuming floating-point arithmetic with a fused multiply-add operation and rounding to nearest, the Cornea-Harrison-Tang method aims to evaluate expressions of the form ab + cd with high relative accuracy. In this article, we provide a rounding error analysis of this method, which unlike previous studies is not restricted to binary floating-point arithmetic but holds for any radix β. We show first that an asymptotically optimal bound on the relative error of this method is 2β u + 2 u 2 /β - 2 u 2 = 2 u + 2/β u 2 + O ( u 3 ), where u = 1/2β 1- p is the unit roundoff in radix β and precision p . Then we show that the possibility of removing the O ( u 2 ) term from this bound is governed by the radix parity and the tie-breaking strategy used for rounding: if β is odd or rounding is to nearest even , then the simpler bound 2 u is obtained, while if β is even and rounding is to nearest away , then there exist floating-point inputs a , b , c , d that lead to a relative error larger than 2 u + 2/β u 2 — 4 u 3 . All these results hold provided underflows and overflows do not occur and under some mild assumptions on p satisfied by IEEE 754-2008 formats.

The interpolation step in the Guruswami-Sudan algorithm is a bivariate interpolation problem with multiplicities commonly solved in the literature using either structured linear algebra or basis reduction of polynomial lattices. This problem has been extended to three or more variables; for this generalization, all fast algorithms proposed so far rely on the lattice approach. In this paper, we reduce this multivariate interpolation problem to a problem of simultaneous polynomial approximations, which we solve using fast structured linear algebra. This improves the best known complexity bounds for the interpolation step of the list-decoding of Reed-Solomon codes, Parvaresh-Vardy codes, and folded Reed-Solomon codes. In particular, for Reed-Solomon list-decoding with re-encoding, our approach has complexity O~(ℓω-1m2(n - k)), where ℓ, m, n, and k are the list size, the multiplicity, the number of sample points, and the dimension of the code, and ω is the exponent of linear algebra; this accelerates the previously fastest known algorithm by a factor of ℓ/m.

This paper deals with the accuracy of complex division in radix-two floating-point arithmetic. Assuming that a fused multiply-add (FMA) instruction is available and that no underflow/overflow occurs, we study how to ensure high relative accuracy in the component wise sense. Since this essentially reduces to evaluating accurately three expressions of the form ac+bd, an obvious approach would be to perform three calls to Kahan's compensated algorithm for 2 by 2 determinants. However, in the context of complex division, two of those expressions are such that ac and bd have the same sign, suggesting that cheaper schemes should be used here (since cancellation cannot occur). We first give a detailed accuracy analysis of such schemes for the sum of two nonnegative products, providing not only sharp bounds on both their absolute and relative errors, but also sufficient conditions for the output of one of them to coincide with the output of Kahan's algorithm. By combining Kahan's algorithm with this particular scheme, we then deduce two new division algorithms. Our first algorithm is a straight-line program whose component wise relative error is always at most 5u+13u2with u the unit round off, we also provide examples of inputs for which the error of this algorithm approaches 5u, thus showing that our upper bound is essentially the best possible. When tests are allowed we show with a second algorithm that the bound above can be further reduced to 4.5u+9u2, and that this improved bound is reasonably sharp.

Graphics and signal processing applications often require that sines and cosines be evaluated at a same floating-point argument, and in such cases a very fast computation of the pair of values is desirable. This paper studies how 32-bit VLIW integer architectures can be exploited in order to perform this task accurately for IEEE single precision (including subnormals). We describe software implementations for sinf, cosf, and sincosf over [-pi/4, pi/4]that have a proven 1-ulp accuracy and whose latency on STMicroelectronics' ST231 VLIW integer processor is 19, 18, and 19 cycles, respectively. Such performances are obtained by introducing a novel algorithm for simultaneous sine and cosine that combines univariate and bivariate polynomial evaluation schemes.

We consider the problem of computing IEEE floating-point squares by means of integer arithmetic. We show how to exploit the specific properties of squaring in order to design and implement algorithms that have much lower latency than those for general multiplication, while still guaranteeing correct rounding. Our algorithms are parameterized by the floating-point format, aim at high instruction-level parallelism (ILP) exposure, and cover all rounding modes. We show further that their C implementation for the binary32 format yields efficient codes for targets like the ST231 VLIW integer processor from ST Microelectronics, with a latency at least 1.75x smaller than that of general multiplication in the same context.

In this paper, we show how to reduce the computation of correctly rounded square roots of binary floating-point data to the fixed-point evaluation of some particular integer polynomials in two variables. By designing parallel and accurate evaluation schemes for such bivariate polynomials, we show further that this approach allows for high instruction-level parallelism (ILP) exposure, and thus, potentially low-latency implementations. Then, as an illustration, we detail a C implementation of our method in the case of IEEE 754-2008 binary32 floating-point data (formerly called single precision in the 1985 version of the IEEE 754 standard). This software implementation, which assumes 32-bit unsigned integer arithmetic only, is almost complete in the sense that it supports special operands, subnormal numbers, and all rounding-direction attributes, but not exception handling (that is, status flags are not set). Finally, we have carried out experiments with this implementation on the ST231, an integer processor from the STMicroelectronics' ST200 family, using the ST200 family VLIW compiler. The results obtained demonstrate the practical interest of our approach in that context: for all rounding-direction attributes, the generated assembly code is optimally scheduled and has indeed low latency (23 cycles).

When implementing a function f in floating-point arithmetic, if we wish correct rounding and good performance, it is important to know if there are input floating-point values x such that f(x) is either the middle of two consecutive floating-point numbers (assuming rounded-to-nearest arithmetic), or a floating-point number (assuming rounded toward ± ∞ or toward 0 arithmetic). In the first case, we say that f(x) is a midpoint, and in the second case, we say that f(x) is an exact point. For some usual algebraic functions and various floating-point formats, we prove whether or not there exist midpoints or exact points. When there exist midpoints or exact points, we characterize them or list all of them (if there are not too many). The results and the techniques presented in this paper can be used in particular to deal with both the binary and the decimal formats defined in the IEEE 754-2008 standard for floating-point arithmetic.

The asymptotically fastest known divide-and-conquer methods for inverting dense structured matrices are essentially variations or extensions of the Morf/Bitmead-Anderson algorithm. Most of them must deal with the growth in length of intermediate generators, and this is done by incorporating various generator compression techniques into the algorithms. One exception is an algorithm by Cardinal, which in the particular case of Cauchy-like matrices avoids such growth by focusing on well-specified, already compressed generators of the inverse. In this paper, we extend Cardinal's method to a broader class of structured matrices including those of Vandermonde, Hankel, and Toeplitz types. Besides, some first experimental results illustrate the practical interest of the approach.

This paper deals with the design and implementation of low latency software for binary floating-point division with correct rounding to nearest. The approach we present here targets a VLIW integer processor of the ST200 family, and is based on fast and accurate programs for evaluating some particular bivariate polynomials. We start by giving approximation and evaluation error conditions that are sufficient to ensure correct rounding. Then we describe the heuristics used to generate such evaluation programs, as well as those used to automatically validate their accuracy. Finally, we propose, for the binary32 format, a complete C implementation of the resulting division algorithm. With the ST200 compiler and compared to previous implementations, the speed-up observed with our approach is by a factor of almost 1.8.

Linear systems with structures such as Toeplitz-, Vandermonde-or Cauchy-likeness can be solved in O~(α2n) operations, where n is the matrix size, α is its displacement rank, and O~denotes the omission of logarithmic factors. We show that for Toeplitz-like and Vandermonde-like trices, this cost can be reduced to O~(αω--1 n), where ω is a feasible exponent for matrix multiplication over the base field. The best known estimate for ω is ω< 2.38, resulting in costs of order O~(α1.38n). We also present consequences for Hermite-Padé approximation and bivariate interpolation.

We study the link between the complexity of polynomial matrix multiplication and the complexity of solving other basic linear algebra problems on polynomial matrices. By polynomial matrices we mean ntimes n matrices in K[x] of degree bounded by d, with K a commutative field. Under the straight-line program model we show that multiplication is reducible to the problem of computing the coefficient of degree d of the determinant. Conversely, we propose algorithms for minimal approximant computation and column reduction that are based on polynomial matrix multiplication; for the determinant, the straight-line program we give also relies on matrix product over K[x] and provides an alternative to the determinant algorithm of [16, 17]. We further show that all these problems can be solved in particular in O (ω) operations in K. Here the "soft O" notation O indicates some missing log (nd) factors and ω is the exponent of matrix multiplication over K.

We show that every perturbation A(λ, ε) of an n x n matrix polynomial A(λ) such that det A(λ) = λm with m ≤ n can be reduced by equivalence transforms to a perturbed matrix polynomial whose leading matrix has maximal Smith form. This yields a reduced form for square perturbed matrix polynomials from which one can easily recover all the eigenvalue leading terms of the form μεβ with β-1 ∈ ℕ*.

In this article, we present an algorithmic approach to the eigenvalue perturbation problem. We show that any matrix perturbation A(ε) of an arbitrary nilpotent Jordan canonical form J with all eigenvalues having an order of the form O(ε1/(a positive integer)) is similar to a matrix perturbation Atilde;(ε) in Arnold normal form that can be seen as generic. Calling A(ε) a κ-matrix and Atilde;(ε) a Lidskii-Arnold matrix, we also provide a reduction algorithm for the computation of the Lidskii-Arnold form of a κ-matrix. It is based on the minimization of the leading Jordan structure J and on Lidskii's genericity conditions for perturbed eigenvalues.

Article A reduction algorithm for matrices depending on a parameter Share on Authors: C.-P. Jeannerod LMC-IMAG, 51 Rue des Mathématiques, 38041 Grenoble Cedex 9, France LMC-IMAG, 51 Rue des Mathématiques, 38041 Grenoble Cedex 9, FranceView Profile , E. Pflügel LMC-IMAG, 51 Rue des Mathématiques, 38041 Grenoble Cedex 9, France LMC-IMAG, 51 Rue des Mathématiques, 38041 Grenoble Cedex 9, FranceView Profile Authors Info & Claims ISSAC '99: Proceedings of the 1999 international symposium on Symbolic and algebraic computationJuly 1999 Pages 121–128https://doi.org/10.1145/309831.309884Online:01 July 1999Publication History 5citation231DownloadsMetricsTotal Citations5Total Downloads231Last 12 Months2Last 6 weeks1 Get Citation AlertsNew Citation Alert added!This alert has been successfully added and will be sent to:You will be notified whenever a record that you have chosen has been cited.To manage your alert preferences, click on the button below.Manage my AlertsNew Citation Alert!Please log in to your account Save to BinderSave to BinderCreate a New BinderNameCancelCreateExport CitationPublisher SiteGet Access