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).

Assume we use a binary floating-point arithmetic and that RN is the round-to-nearest function. Also assume that c is a constant or a real function of one or more variables, and that we have at our disposal a correctly rounded implementation of c, say ĉ = RN(c). For evaluating x • c (resp. x/c or c/x), the natural way is to replace it by RN(x • ĉ) (resp. RN(x/ĉ) or RN(ĉ/x)), that is, to call function ĉ and to perform a floating-point multiplication or division. This can be generalized to the approximation of n/d by RN(n̂ / d̂ ) and the approximation of n • d by RN(n̂ • d̂ ), where n̂ = RN(n) and d̂ = RN(d), and n and d are functions for which we have at our disposal a correctly rounded implementation. We discuss tight error bounds in ulps of such approximations. From our results, one immediately obtains tight error bounds for calculations such as x * pi, ln(2)/x, x/(y + z), (x + y) * z, x/sqrt(y), sqrt(x)/y, (x + y)(z + t), (x + y)/(z + t), (x + y)/(zt), etc. in floating-point arithmetic.

The computation of Fast Fourier Transforms (FFTs) in floating-point arithmetic is inexact due to roundings, and for some applications it can prove very useful to know a tight bound on the final error. Although it can be almost attained by specifically built input values, the best known error bound for the Cooley-Tukey FFT seems to be much larger than most actually obtained errors. Also, interval arithmetic can be used to compute a bound on the error committed with a given set of input values, but it is in general considered hampered with large overestimation. We report results of intensive computations to test the two approaches, in order to estimate the numerical performance of state-of-the-art bounds. Surprisingly enough, we observe that while interval arithmetic-based bounds are overestimated, they remain, in our computations, tighter than general known bounds.

We consider the computation of the Euclidean (or L2) norm of an n -dimensional vector in floating-point arithmetic. We review the classical solutions used to avoid spurious overflow or underflow and/or to obtain very accurate results. We modify a recently published algorithm (that uses double-word arithmetic) to allow for a very accurate solution, free of spurious overflows and underflows. To that purpose, we use a double-word square-root algorithm of which we provide a tight error analysis. The returned L2 norm will be within very slightly more than 0.5 ulp from the exact result, which means that we will almost always provide correct rounding.

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.

Recently, a complete set of algorithms for manipulating double-word numbers (some classical, some new) was analyzed [ 16 ]. We have formally proven all the theorems given in that article, using the Coq proof assistant. The formal proof work led us to: (i) locate mistakes in some of the original paper proofs (mistakes that, however, do not hinder the validity of the algorithms), (ii) significantly improve some error bounds, and (iii) generalize some results by showing that they are still valid if we slightly change the rounding mode. The consequence is that the algorithms presented in [ 16 ] can be used with high confidence, and that some of them are even more accurate than what was believed before. This illustrates what formal proof can bring to computer arithmetic: beyond mere (yet extremely useful) verification, correction, and consolidation of already known results, it can help to find new properties. All our formal proofs are freely available.

Expressions such as$ax^{2}, axy$, or$ax^{3}$, where$a$is a constant, are not unfrequent in computing. There are several ways of parenthesizing them (and therefore, choosing the order of evaluation). Depending on the value of$a$, is there a more accurate evaluation order? We discuss this point (with a small digression on spurious underflows and overflows).

The 2019 version of the IEEE 754 Standard for Floating-Point Arithmetic recommends that new “augmented” operations should be provided for the binary formats. These operations use a new “rounding direction”: round-to-nearestties-to-zero. We show how they can be implemented using the currently available operations, using round-to-nearestties-to-evenwith a partial formal proof of correctness.

We introduce algorithms for splitting a positive binary floating-point number into two numbers of around half the system precision, using arithmetic operations all rounded either toward -∞ or toward +∞. We use these algorithms to compute “exact” products (i.e., to express the product of two floating-point numbers as the unevaluated sum of two floating-point numbers, the rounded product and an error term). This is similar to the classical Dekker product, adapted here to directed roundings.

Quaternions form a set of four global but not unique parameters, which can represent three-dimensional rotations in a non-singular way. They are frequently used in computer graphics, drone and aerospace vehicle control. Floating-point quaternion operations (addition, multiplication, reciprocal, norm) are often implemented “by the book”. Although all usual implementations are algebraically equivalent, their numerical behavior can be quite different. For instance, the arithmetic operations on quaternions as well as conversion algorithms to/from rotation matrices are subject to spurious under/overflow (an intermediate calculation underflows or overflows, making the computed final result irrelevant, although the exact result is in the domain of the representable numbers). The goal of this paper is to analyze and then propose workarounds and better accuracy alternatives for such algorithms.

In this article, we review some of the classical methods used for quickly obtaining low-precision approximations to the elementary functions. Then, for each of the three main classes of elementary function algorithms (shift-and-add algorithms, polynomial or rational approximations, and table-based methods) and for the additional, specific to approximate computing, “bit-manipulation” techniques, we examine what can be done for obtaining very fast estimates of a function, at the cost of a (controlled) loss in terms of accuracy.

We are interested in obtaining error bounds for the classical Cooley-Tukey fast Fourier transform algorithm in floating-point arithmetic, for the 2-norm as well as for the infinity norm. For that purpose, we also give some results on the relative error of the complex multiplication by a root of unity, and on the largest value that can take the real or imaginary part of one term of the fast Fourier transform of a vector x , assuming that all terms of x have real and imaginary parts less than some value b .

We deal with accurate complex multiplication in binary floating-point arithmetic, with an emphasis on the case where one of the operands is a "double-word" number. We provide an algorithm that returns a complex product with normwise relative error bound close to the best possible one, i.e., the rounding unit u.

The normal and complementary error functions are ubiquitous special functions for any mathematical library. They have a wide range of applications. Practical applications call for customized implementations that have strict accuracy requirements. Accurate numerical implementation of these functions is, however, non-trivial. In particular, the complementary error function erfc for large positive arguments heavily suffers from cancellation, which is largely due to its asymptotic behavior. We provide a semi-automatic code generator for the erfc function which is parameterized by the user-given bound on the relative error. Our solution exploits the asymptotic expression of erfc and leverages the automatic code generator Metalibm that provides accurate polynomial approximations. A fine-grained a priori error analysis provides a libm developer with the required accuracy for each step of the evaluation. In critical parts, we exploit double-word arithmetic to achieve implementations that are fast, yet accurate up to 50 bits, even for large input arguments. We demonstrate that for high required accuracies the automatically generated code has performance comparable to that of the standard libm and for lower ones our code demonstrated roughly 25% speedup.

Triple-word arithmetic consists in representing high-precision numbers as the unevaluated sum of three floating-point numbers (with “nonoverlapping” constraints that are explicited in the paper). We introduce and analyze various algorithms for manipulating triple-word numbers: rounding a triple-word number to a floating-point number, adding, multiplying, dividing, and computing square-roots of triple-word numbers, etc. We compare our algorithms, implemented in the Campary library, with other solutions of comparable accuracy. It turns out that our new algorithms are significantly faster than what one would obtain by just using the usual floating-point expansion algorithms in the special case of expansions of length 3.

We present an automatic method for the evaluation of functions via polynomial or rational approximations and its hardware implementation, on FPGAs. These approximations are evaluated using Ercegovac's iterative E-method adapted for FPGA implementation. The polynomial and rational function coefficients are optimized such that they satisfy the constraints of the E-method. We present several examples of practical interest; in each case a resource-efficient approximation is proposed and comparisons are made with alternative approaches.

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.

Semidefinite programming (SDP) is widely used in optimization problems with many applications, however, certain SDP instances are ill-posed and need more precision than the standard double-precision available. Moreover, these problems are large-scale and could benefit from parallelization on specialized architectures such as GPUs. In this article, we implement and evaluate the performance of a floating-point expansion-based arithmetic library (CAMPARY) in the context of such numerically highly accurate SDP solvers. We plugged-in CAMPARY with the state-of-the-art SDPA solver for both CPU and GPU-tuned implementations. We compare and contrast both the numerical accuracy and performance of SDPA-GMP, -QD and -DD, which employ other multiple-precision arithmetic libraries against SDPA-CAMPARY. We show that CAMPARY is a very good trade-off for accuracy and speed when solving ill-conditioned SDP problems.

The papers in this special issue focus on computer arithmetic which is used in many applications, usually totally silently (one should keep in mind that even when running programs that are not at all numeric, memory addresses are computed, which involves additions, multiplications, and sometimes divisions). However, in some areas, it plays a central role.

The 2Sum and Fast2Sum algorithms are important building blocks in numerical computing. They are used (implicitely or explicitely) in many compensated algorithms (such as compensated summation or compensated polynomial evaluation). They are also used for manipulating floating-point expansions . We show that these algorithms are much more robust than it is usually believed: The returned result makes sense even when the rounding function is not round-to-nearest, and they are almost immune to overflow.

We analyze several classical basic building blocks of double-word arithmetic (frequently called “double-double arithmetic” in the literature): the addition of a double-word number and a floating-point number, the addition of two double-word numbers, the multiplication of a double-word number by a floating-point number, the multiplication of two double-word numbers, the division of a double-word number by a floating-point number, and the division of two double-word numbers. For multiplication and division we get better relative error bounds than the ones previously published. For addition of two double-word numbers, we show that the previously published bound was incorrect, and we provide a new relative error bound. We introduce new algorithms for division. We also give examples that illustrate the tightness of our bounds.

Elementary functions from the mathematical library input and output floating-point numbers. However it is possible to implement them purely using integer/fixed-point arithmetic. This option was not attractive between 1985 and 2005, because mainstream processor hardware supported 64-bit floating-point, but only 32-bit integers. This has changed in recent years, in particular with the generalization of native 64-bit integer support. The purpose of this article is therefore to reevaluate the relevance of computing floating-point functions in fixed-point. For this, several variants of the double-precision logarithm function are implemented and evaluated. Formulating the problem as a fixed-point one is easy after the range has been (classically) reduced. Then, 64-bit integers provide slightly more accuracy than 53-bit mantissa, which helps speed up the evaluation. Finally, multi-word arithmetic, critical for accurate implementations, is much faster in fixed-point, and natively supported by recent compilers. Thanks to all this, a purely integer implementation of the correctly rounded double-precision logarithm outperforms the previous state of the art, with the worst-case execution time reduced by a factor 5. This work also introduces variants of the logarithm that input a floating-point number and output the result in fixed-point. These are shown to be both more accurate and more efficient than the traditional floating-point functions for some applications.

Some important computational problems must use a floating-point (FP) precision several times higher than the hardware-implemented available one. These computations critically rely on software libraries for high-precision FP arithmetic. The representation of a high-precision data type crucially influences the corresponding arithmetic algorithms. Recent work showed that algorithms for FP expansions, that is, a representation based on unevaluated sum of standard FP types, benefit from various high-performance support for native FP, such as low latency, high throughput, vectorization, threading, etc. Bailey's QD library and its corresponding Graphics Processing Unit (GPU) version, GQD, are such examples. Despite using native FP arithmetic as the key operations, QD and GQD algorithms are focused on double-double or quad-double representations and do not generalize efficiently or naturally to a flexible number of components in the FP expansion. In this paper, we introduce a new multiplication algorithm for FP expansion with flexible precision, up to the order of tens of FP elements in mind. The main feature consists in the partial products being accumulated in a special designed data structure that has the regularity of a fixed-point representation while allowing the computation to be naturally carried out using native FP types. This allows us to easily avoid unnecessary computation and to present rigorous accuracy analysis transparently. The algorithm, its correctness and accuracy proofs and some performance comparisons with existing libraries are all contributions of this paper.

GPUs are an important hardware development platform for problems where massive parallel computations are needed. Many of these problems require a higher precision than the standard double floating-point (FP) available. One common way of extending the precision is the multiple-component approach, in which real numbers are represented as the unevaluated sum of several standard machine precision FP numbers. This representation is called a FP expansion and it offers the simplicity of using directly available and highly optimized FP operations. In this article we present new data-parallel algorithms for adding and multiplying FP expansions specially designed for extended precision computations on GPUs. These are generalized algorithms that can manipulate FP expansions of different sizes (from double-double up to a few tens of doubles) and ensure a certain worst case error bound on the results.

We introduce an algorithm to compare a binary floating-point (FP) number and a decimal FP number, assuming the “binary encoding” of the decimal formats is used, and with a special emphasis on the basic interchange formats specified by the IEEE 754-2008 standard for FP arithmetic. It is a two-step algorithm: a first pass, based on the exponents only, quickly eliminates most cases, then, when the first pass does not suffice, a more accurate second pass is performed. We provide an implementation of several variants of our algorithm, and compare them.

Many numerical problems require a higher computing precision than the one offered by standard floating-point (FP) formats. One common way of extending the precision is to represent numbers in amultiple componentformat. By using the so-calledfloating-point expansions, real numbers are represented as the unevaluated sum of standard machine precision FP numbers. This representation offers the simplicity of using directly available, hardware implemented and highly optimized, FP operations. It is used by multiple-precision libraries such as Bailey's QD or the analogue Graphics Processing Units (GPU) tuned version, GQD. In this article we briefly revisit algorithms for adding and multiplying FP expansions, then we introduce and prove new algorithms for normalizing, dividing and square rooting of FP expansions. The new method used for computing the reciprocal${a}^{-1}$and the square root$\sqrt{a}$of a FP expansion$a$is based on an adapted Newton-Raphson iteration where the intermediate calculations are done using “truncated” operations (additions, multiplications) involving FP expansions. We give here a thorough error analysis showing that it allows very accurate computations. More precisely, after$q$iterations, the computed FP expansion$x=x_0+\ldots +x_{2^q-1}$satisfies, for the reciprocal algorithm, the relative error bound:$\left|({x-a^{-1}})/{a^{-1}}\right| \le 2^{-2^q(p-3)-1}$and, respectively, for the square root one:$\left|x-{1}/{\sqrt{a}}\right| \le {2^{-2^q(p-3)-1}}/{\sqrt{a}}$, where$p> 2$is the precision of the FP representation used ($p=24$for single precision and$p=53$for double precision).

In their book, Scientific Computing on the Itanium , Cornea et al. [2002] introduce an accurate algorithm for evaluating expressions of the form ab + cd in binary floating-point arithmetic, assuming an FMA instruction is available. They show that if p is the precision of the floating-point format and if u = 2 -p , the relative error of the result is of order u . We improve their proof to show that the relative error is bounded by 2 u +7 u 2 +6 u 3 . Furthermore, by building an example for which the relative error is asymptotically (as p → ∞ or, equivalently, as u → 0) equivalent to 2 u , we show that our error bound is asymptotically optimal.

Many numerical problems require a higher computing precision than that offered by common floating point (FP) formats. One common way of extending the precision is to represent numbers in a multiple component format. With so-called floating point expansions, numbers are represented as the unevaluated sum of standard machine precision FP numbers. This format offers the simplicity of using directly available and highly optimized FP operations and is used by multiple-precisions libraries such as Bailey's QD or the analogue Graphics Processing Units tuned version, GQD. In this article we present a new algorithm for computing the reciprocal FP expansion a-1of a FP expansion a. Our algorithm is based on an adapted Newton-Raphson iteration where we use “truncated” operations (additions, multiplications) involving FP expansions. The thorough error analysis given shows that our algorithm allows for computations of very accurate quotients. Precisely, after q ≤ 0 iterations, the computed FP expansion x = x0+ ... + x2q-1satisfies the relative error bound |x-a-1/a-1|≤2-2q(p-3)-1, where p > 2 is the precision of the FP representation used (p = 24 for single precision and p = 53 for double precision).

We introduce a software-oriented algorithm that allows one to quickly compare a binary64 floating-point (FP) number and a decimal64 FP number, assuming the "binary encoding" of the decimal formats specified by the IEEE 754-2008 standard for FP arithmetic is used. It is a two-step algorithm: a first pass, based on the exponents only, makes it possible to quickly eliminate most cases, then when the first pass does not suffice, a more accurate second pass is required. We provide an implementation of several variants of our algorithm, and compare them.

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.

A very simple test, introduced by Ziv, allows one to determine if an approximation to the value f(x) of an elementary function at a given point x suffices to return the floating-point number nearest f(x) . The same test may be used when implementing floating-point operations with input and output operands of different formats, using arithmetic operators tailored for manipulating operands of the same format. That test depends on a “magic constant” e . We show how to choose that constant e to make the test reliable and efficient. Various cases are considered, depending on the availability of an fma instruction, and on the range of f(x) .

We present a new way of approximating the sine and cosine functions by a few table look-ups and additions. It consists in first reducing the input range to a very small interval by using rotations with "(M, p, k) friendly angles", proposed in this work, and then by using a bipartite table method ina small interval. An implementation of the method for 24-bit case is described and compared with CORDIC. Roughly, the proposed scheme offers a speedup of 2 compared with an unfolded double-rotation radix-2 CORDIC.

This paper presents a study of some basic blocks needed in the design of floating-point summation algorithms. In particular, in radix-2 floating-point arithmetic, we show that among the set of the algorithms with no comparisons performing only floating-point additions/subtractions, the 2Sum algorithm introduced by Knuth is minimal, both in terms of number of operations and depth of the dependency graph. We investigate the possible use of another algorithm, Dekker's Fast2Sum algorithm, in radix-10 arithmetic. We give methods for computing, in radix 10, the floating-point number nearest the average value of two floating-point numbers. We also prove that under reasonable conditions, an algorithm performing only round-to-nearest additions/subtractions cannot compute the round-to-nearest sum of at least three floating-point numbers. Starting from an algorithm due to Boldo and Melquiond, we also present new results about the computation of the correctly-rounded sum of three floating-point numbers. For a few of our algorithms, we assume new operations defined by the recent IEEE 754-2008 Standard are available.

Define an "augmented precision" algorithm as an algorithm that returns, in precision-p floating-point arithmetic, its result as the unevaluated sum of two floating-point numbers, with a relative error of the order of 2-2p. Assuming an FMA instruction is available, we perform a tight error analysis of an augmented precision algorithm for the square root, and introduce two slightly different augmented precision algorithms for the 2D-norm √x2+y2. Then we give tight lower bounds on the minimum distance (in ulps) between √x2+y2and a midpoint when √x2+y2is not itself a midpoint. This allows us to determine cases when our algorithms make it possible to return correctly-rounded 2D-norms.

Solving the Table Maker's Dilemma, for a given function and a given target floating-point format, requires testing the value of the function, with high precision, at a very large number of consecutive values. We give an algorithm that allows for performing such computations on a very regular architecture, and present an FPGA implementation of that algorithm.

The fused multiply accumulate-add (FMA) instruction, specified by the IEEE 754-2008 Standard for Floating-Point Arithmetic, eases some calculations, and is already available on some current processors such as the Power PC or the Itanium. We first extend an earlier work on the computation of the exact error of an FMA (by giving more general conditions and providing a formal proof). Then, we present a new algorithm that computes an approximation to the error of an FMA, and provide error bounds and a formal proof for that algorithm.

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.

During any composite computation, there is a constant need for rounding intermediate results before they can participate in further processing. Recently, a class of number representations denoted RN-Codings were introduced, allowing an unbiased rounding-to-nearest to take place by a simple truncation, with the property that problems with double-roundings are avoided. In this paper, we first investigate a particular encoding of the binary representation. This encoding is generalized to any radix and digit set; however, radix complement representations for even values of the radix turn out to be particularly feasible. The encoding is essentially an ordinary radix complement representation with an appended round-bit, but still allowing rounding-to-nearest by truncation, and thus avoiding problems with double-roundings. Conversions from radix complement to these round-to-nearest representations can be performed in constant time, whereas conversion the other way, in general, takes at least logarithmic time. Not only is rounding-to-nearest a constant time operation, but so is also sign inversion, both of which are at best log-time operations on ordinary two's complement representations. Addition and multiplication on such fixed-point representations are first analyzed and defined in such a way that rounding information can be carried along in a meaningful way, at minimal cost. The analysis is carried through for a compact (canonical) encoding using two's complement representation, supplied with a round-bit. Based on the fixed-point encoding, it is shown possible to define floating-point representations, and a sketch of the implementation of an FPU is presented.

We propose algorithms and provide some related results that make it possible to implement decimal floating-point arithmetic on a processor that does not have decimal operators, using the available binary floating-point functions. In this preliminary study, we focus on round-to-nearest mode only. We show that several functions in decimal32 and dec-imal64 arithmetic can be implemented using binary64 and binaryl28 floating-point arithmetic, respectively. We discuss the decimal square root and some transcendental functions. We also consider radix conversion algorithms.

Since the introduction of the Fused Multiply and Add (FMA) in the IEEE-754-2008 standard for floatingpoint arithmetic, division based on Newton-Raphson's iterations becomes a viable alternative to SRT-based divisions. The Newton-Raphson iterations were already used in some architecture prior to the revision of the IEEE-754 norm. For example, Itanium architecture already used this kind of iterations. Unfortunately, the proofs of the correctness of binary algorithms do not extend to the case of decimal floating-point arithmetic. In this paper, we present general methods to prove the correct rounding of division algorithms using Newton-Raphson's iterations in software, for radix 2 and radix 10 floating-point arithmetic.

We introduce several algorithms for accurately evaluating powers to a positive integer in floating-point arithmetic, assuming a fused multiply-add (fma) instruction is available. For bounded, yet very large values of the exponent, we aim at obtaining correctly rounded results in round-to-nearest mode, that is, our algorithms return the floating-point number that is nearest the exact value.

This paper presents a study of some basic blocks needed in the design of floating-point summation algorithms. In particular, we show that among the set of the algorithms with no comparisons performing only floating-point additions/subtractions, the 2Sum algorithm introduced by Knuth is minimal, both in terms of number of operations and depth of the dependency graph. Under reasonable conditions, we also prove that no algorithms performing only round-to-nearest additions/subtractions exist to compute the round-to-nearest sum of at least three floating-point numbers. Starting from an algorithm due to Boldo and Melquiond, we also present new results about the computation of the correctly-rounded sum of three floating-point numbers.

We present a design and implementation of a radix-4 complex division unit with prescaling of the operands. Specifically, we extend the treatment of the residual bound and errors due to the use of truncated redundant representation. The requirements for prescaling tables are simplified and a detailed specification of the table design is given. All principal components used in the design are described and the proposed optimizations are explained. The target platform for implementation was an Altera Stratix II FPGA for which we report timing and area requirements. For a precision of 36 bits, the implementation uses 1185 ALUTs, achieving a latency of 157 ns. The maximum clock frequency is 173.49 MHz.

In this paper we extend the domain of applicability of the E-method [7, 8], as a hardware-oriented method for evaluating elementary functions using polynomial and rational function approximations. The polynomials and rational functions are computed by solving a system of linear equations using digit-serial iterations on simple and highly regular hardware. For convergence, these systems must be diagonally dominant. The E-method offers an efficient way for the fixed-point evaluation of polynomials and rational functions if their coefficients conform to the diagonal dominance condition. Until now, there was no systematic approach to obtain good approximations to f over an interval [a, b] by rational functions satisfying the constraints required by the E-method. In this paper, we present such an approach which is based on linear programming and lattice basis reduction. We also discuss a design and performance characteristics of a corresponding implementation.

Reconfigurable circuits now have a capacity that allows them to be used as floating-point accelerators. They offer massive parallelism, but also the opportunity to design optimised floating-point hardware operators not available in microprocessors. Multiplication by a constant is an important example of such an operator. This article presents an architecture generator for the correctly rounded multiplication of a floating-point number by a constant. This constant can be a floating-point value, but also an arbitrary irrational number. The multiplication of the significands is an instance of the well-studied problem of constant integer multiplication, for which improvement to existing algorithms are also proposed and evaluated.

Since redundant number systems allow constant time addition, they are often at the heart of modular multipliers designed for public key cryptography (PKC) applications. Indeed, PKC involves large operands (160 to 1024 bits) and several researchers proposed carry-save or borrow-save algorithms. However, these number systems do not take advantage of the dedicated carry logic available in modern field programmable gate arrays (FPGAs). To overcome this problem, we suggest to perform modular multiplication in a high-radix carry-save number system, where a sum bit of the carry-save representation is replaced by a sum word. Two digits are then added by means of a small Carry-Ripple Adder (CRA). Furthermore, we propose an algorithm which selects the best high-radix carry-save representation for a given modulus, and generates a synthesizable VHDL description of the operator.

We introduce an algorithm for multiplying a floating-point number x by a constant C that is not exactly representable in floating-point arithmetic. Our algorithm uses a multiplication and a fused multiply and add instruction. Such instructions are available in some modern processors such as the IBM Power PC and the Intel/HP Itanium. We give three methods for checking whether, for a given value of C and a given floating-point format, our algorithm returns a correctly rounded result for any x. When it does not, some of our methods return all of the values x for which the algorithm fails. The three methods are complementary: The first two do not always allow one to conclude, yet they are simple enough to be used at compile time, while the third one always either proves that our algorithm returns a correctly rounded result for any x or gives all of the counterexamples. We generalize our study to the case where a wider internal format is used for the intermediate calculations, which gives a fourth method. Our programs and some additional information (such as the case where an arbitrary nonbinary even radix is used), as well as examples of runs of our programs, can be downloaded from http://perso.ens-lyon.fr/iean-michel.muller/MultConstant.html.

A hardware-oriented method for evaluating complex polynomials by solving iteratively a system of linear equations is proposed. Its implementation uses a digit-serial iterations on simple and highly regular hardware. The operations involved are defined over the reals. We describe a complex-to-real transform, a complex polynomial evaluation algorithm, the convergence conditions, and a corresponding design and implementation. The latency and the area are estimated for the radix-2 case. The main features of the method are: the latency of about m cycles for an m-bit precision; the cycle time independent of the precision; a design consisting of identical modules; and a digit-serial connections between the modules. The number of modules, each roughly corresponding to serial-parallel multiplier without a carry-propagate adder, is 2(n + I) for evaluating an n-th degree complex polynomial. The method can also be used to compute all successive integer powers of the complex argument with the same latency and a similar implementation cost.

Redundant number representations are generally used to allow constant time additions, based on the fact that only bounded carry-ripples take place. But, carries may ripple out into positions which may not be needed to represent the final value of the result and, thus, a certain amount of leading guard digits are needed to correctly determine the result. Also, when cancellation during subtractions occurs, there may be nonzero digits in positions not needed to represent the result of the calculation. It is shown here that, for normal redundant digit sets with radix greater than two, a single guard digit is sufficient to determine the value of such an arbitrary length prefix of leading nonzero digits. This is also the case for the unsigned carry-save representation, whereas two guard digits are sufficient, and may be necessary, for additions in the binary signed-digit and 2's complement carry-save representations. Thus, only the guard digits need to be retained during sequences of additions and subtractions. At suitable points, the guard digits may then be converted into a single digit, representing the complete prefix.

Polynomial approximations are almost always used when implementing functions on a computing system. In most cases, the polynomial that best approximates (for a given distance and in a given interval) a function has coefficients that are not exactly representable with a finite number of bits. And yet, the polynomial approximations that are actually implemented do have coefficients that are represented with a finite---and sometimes small---number of bits. This is due to the finiteness of the floating-point representations (for software implementations), and to the need to have small, hence fast and/or inexpensive, multipliers (for hardware implementations). We then have to consider polynomial approximations for which the degree- i coefficient has at most m i fractional bits; in other words, it is a rational number with denominator 2 m i . We provide a general and efficient method for finding the best polynomial approximation under this constraint. Moreover, our method also applies if some other constraints (such as requiring some coefficients to be equal to some predefined constants or minimizing relative error instead of absolute error) are required.

The fused multiply accumulate instruction (fused-mac) that is available on some current processors such as the Power PC or the Itanium eases some calculations. We give examples of some floating-point functions (such as ulp(x) or Nextafter(x, y)), or some useful tests, that are easily computable using a fused-mac. Then, we show that, with rounding to the nearest, the error of a fused-mac instruction is exactly representable as the sum of two floating-point numbers. We give an algorithm that computes that error.

We introduce an algorithm for multiplying a floating-point number x by a constant C that is not exactly representable in floating-point arithmetic. Our algorithm uses a multiplication and a fused multiply and add instruction. We give methods for checking whether, for a given value of C and a given floating-point format, our algorithm returns a correctly rounded result for any x. When it does not, our methods give the values x for which it does not.

Algorithms for Euclidean (i.e., integer) division by a constant operation are presented. They allow fast computation for some values of the divisor (known at compile time) or also when both quotient and modulus are required. These algorithms are based on the multiply-accumulate instruction and the 40-bit arithmetic available in DSPs such as the ST100 DSP from STMicroelectronics. The results are demonstrated in the case of standard speech coding applications.

The RN-codings, where "RN" stands for "round to nearest", are particular cases of signed digit representations, for which rounding to nearest is always identical to truncation. In radix 2, booth recoding is an RN-coding. In this paper, we suggest several multiplication algorithms able to handle RN-codings, and we analyze their properties.

We propose a digit recurrence algorithm for division in real and complex number domains using a variable radix. The objective of the approach is to simplify the prescaling of the operands by using a suitable low radix, and switch to higher radices in the remaining iterations to reduce their number. The prescaling is used to allow a simple quotient digit selection by rounding of the residual. We discuss the algorithm, its implementation, and estimate its time and cost characteristics with respect to fixed high radix division algorithms.

Range-reduction is a key point for getting accurate elementary function routines. We introduce a new algorithm that is fast for input arguments belonging to the most common domains, yet accurate over the full double-precision range.

A table-based method for high-speed function approximation in single-precision floating-point format is presented in this paper. Our focus is the approximation of reciprocal, square root, square root reciprocal, exponentials, logarithms, trigonometric functions, powering (with a fixed exponent p), or special functions. The algorithm presented here combines table look-up, an enhanced minimax quadratic approximation, and an efficient evaluation of the second-degree polynomial (using a specialized squaring unit, redundant arithmetic, and multioperand addition). The execution times and area costs of an architecture implementing our method are estimated, showing the achievement of the fast execution times of linear approximation methods and the reduced area requirements of other second-degree interpolation algorithms. Moreover, the use of an enhanced minimax approximation which, through an iterative process, takes into account the effect of rounding the polynomial coefficients to a finite size allows for a further reduction in the size of the look-up tables to be used, making our method very suitable for the implementation of an elementary function generator in state-of-the-art DSPs or graphics processing units (GPUs).

We present techniques for accelerating the floating-point computation of x/y when y is known before x. The proposed algorithms are oriented toward architectures with available fused-mac operations. The goal is to get exactly the same result as with usual division with rounding to nearest. It is known that the advanced computation of 1/y allows performing correctly rounded division in one multiplication plus two fused-macs. We show algorithms that reduce this latency to one multiplication and one fused-mac. This is achieved if a precision of at least n+1 bits is available, where n is the number of mantissa bits in the target format, or if y satisfies some properties that can be easily checked at compile-time. This requires a double-word approximation of 1/y (we also show how to get it). Compilers to accelerate some numerical programs without loss of accuracy can use these techniques.

We aim at evaluating elementary and special functions using small tables and small, rectangular, multipliers. To do that, we show how accurate polynomial approximations whose order-1 coefficients are small in size (a few bits only) can be computed. We compare the obtained results with similar work in the recent literature.

We adapt the radix-r digit-recurrence division algorithm to complex division. By prescaling the operands, we make the selection of quotient digits simple. This leads to a simple hardware implementation, and allows correct rounding of complex quotient. To reduce large prescaling tables required for radices greater than 4, we adapt the bipartite-table method to multiple-operand functions.

This paper presents upper bounds on the number of zeros and of ones after the rounding bit for algebraic functions. These functions include reciprocal, division, square root, and reciprocal square root, which have been considered in previous work. We propose simpler proofs for the previously given bounds and generalize to all algebraic functions. We also determine cases for which the bound is achieved for square root. As is mentioned in the previous work, these bounds are useful for determining the precision required in the computation of approximations in order to be able to perform correct rounding. We consider rounding to nearest, but the results can be easily extended to other rounding modes.

We give the results of a four-year search for the worst cases for correct rounding of the major elementary functions in double precision. These results allow the design of reasonably fast routines that will compute these functions with correct rounding, at least in some interval, for any of the four rounding modes specified by the IEEE-754 standard. They will also allow one to easily test libraries that are claimed to provide correctly rounded functions.

A method for the calculation of faithfully rounded single-precision floating-point powering (X/sup p/) is proposed in this paper. This method employs table look-up and a second-degree minimax approximation, which allows the employment of reduced size tables to store the coefficients from the polynomial approximation. A specialized squaring unit and a fused accumulation tree carry out with the computation of the quadratic polynomial. Both unfolded and pipelined architectures are presented, and the results of a pre-layout synthesis performed using CMOS 0.35 /spl mu/m technology are shown, achieving a 50% area reduction from linear approximation methods, and with improved speed over other second-degree approximation based algorithms. The pipelined architecture has a latency of three cycles and a throughput of one result per cycle.

A FPGA implementation of a method for the calculation of faithfully rounded single-precision floating-point powering (X/sup p/) is presented in this paper. A second-degree minimax polynomial approximation is used, together with the employment of table look-up, a specialized squaring unit and a fused accumulation tree. The FPGA implementation of an architecture with a latency of 3 cycles and a throughput of one result per cycle has been performed using a Xilinx XC4036XL device. The implemented unit has an operation frequency over 33 MHz.

The aim of this paper is to accelerate division, square root, and square root reciprocal computations when the Goldschmidt method is used on a pipelined multiplier. This is done by replacing the last iteration by the addition of a correcting term that can be looked up during the early iterations. We describe several variants of the Goldschmidt algorithm, assuming 4-cycle pipelined multiplier, and discuss obtained number of cycles and error achieved. Extensions to other than 4-cycle multipliers are given. If we call G/sub m/ the Goldschmidt algorithm with m iterations, our variants allow us to reach an accuracy that is between that of G/sub 3/ and that of G/sub 4/, with a number of cycle equal to that of G/sub 3/.

This paper deals with the computation of reciprocals, square roots, inverse square roots, and some elementary functions using small tables, small multipliers, and, for some functions, a final "large" (almost full-length) multiplication. We propose a method, based on argument reduction and series expansion, that allows fast evaluation of these functions in high precision. The strength of this method is that the same scheme allows the computation of all these functions. We estimate the delay, the size/number of tables, and the size/number of multipliers and compare with other related methods.

The Table Maker's Dilemma is the problem of always getting correctly rounded results when computing the elementary functions. After a brief presentation of this problem, we present new developments that have helped us to solve this problem for the double-precision exponential function in a small domain. These new results show that this problem can be solved, at least for the double-precision format, for the most usual functions.

We present a new class of number systems, called Semi-Logarithmic Number Systems, that constitute a family of various compromises between floating-point and logarithmic number systems. This allows trade between the speed of the arithmetic operations and the size of the required tables. We give arithmetic algorithms (addition/subtraction, multiplication, division) for the Semi-Logarithmic Number Systems, and we compare these number systems to the classical floating-point or logarithmic number systems.

The Table Maker's Dilemma is the problem of always getting exactly rounded results when computing the elementary functions. After a brief presentation of this problem, we present new developments that helped us to solve this problem for the double precision exponential function in a small domain. These new results show that this problem can be solved, at least for the double precision format, for the most usual functions.

We propose in this paper a new algorithm and architecture for performing divisions in residue number systems. Our algorithm is suitable for residue number systems with large moduli, with the aim of manipulating very large integers on a parallel computer or a special-purpose architecture. The two basic features of our algorithm are on one hand the use of a high-radix division method, and on the other hand the use of a floating-point arithmetic that should run in parallel with the modular arithmetic.

We present a new class of number systems, called semi-logarithmic number systems, that constitute a family of various compromises between floating-point and logarithmic number systems. We propose arithmetic algorithms for the semi-logarithmic number systems, and we compare these number systems to the classical floating-point or logarithmic number systems.>

A new algorithm for computing the complex logarithm and exponential functions is proposed. This algorithm is based on shift-and-add elementary steps, and it generalizes some algorithms by Briggs and De Lugish (1970), as well as the CORDIC algorithm. It can easily be used to compute the classical real elementary functions (sin, cos, arctan, ln, exp). This algorithm is more suitable for computations in a redundant number system than the CORDIC algorithm, since there is no scaling factor when computing trigonometric functions.>

After a short introduction to on-line computing, we prove that the functions computable in on-line by a finite automaton are piecewise affine functions whose coefficients are rational numbers (i.e., the functions f(x)=ax+b, or f(x,y)=ax+by+c where a, b, and c are rational). A consequence of this study is that multiplication, division and elementary functions of operands of arbitrarily long length cannot be performed using bounded-size operators.>

An algorithm for computing complex logarithms and exponentials is proposed. The algorithm is based on shift-and-add elementary steps, and it generalizes the Cordic algorithm. It can compute the usual real elementary functions. This algorithm is more suitable for computations in a redundant number system than Cordic, since there is no scaling factor for computation of trigonometric functions.>

The parallel and online (i.e., digit serial, most significant digit first) evaluation of polynomials and inverses of polynomials is dealt with. New algorithms and architectures are proposed for such evaluations. A 3-D implementation model is presented.>

The authors present a new modular architecture for the online evaluation of power series. It can be used to quickly compute any function that can be approximated by the first terms of its Taylor expansion (i.e., most math functions). For trigonometric functions, the method matches Cordic-like methods and leads to more regular architectures. The authors also present a VLSI implementation of this architecture.>

After a brief survey of the CORDIC algorithm, some new results that allow fast and easy signed-digit implementation of CORDIC, without modifying the basic iteration step are given. A slight modification would make it possible to use a carry-save representation of numbers, instead of a signed-digit one. The method, called the branching CORDIC method, consists of performing in parallel two classic CORDIC rotations. It gives a constant normalization factor. An online implementation of the algorithm is proposed with an online delay equal to 5 for the sine and cosine functions.>

An extension of the CORDIC (coordinate rotation digital computer) algorithm that makes it possible to compute the functions cos/sup -1/, sin/sup -1/, square root 1-t/sup 2/, sinh/sup -1/ cosh/sup -1/, and square root 1+t/sup 2/ is presented. The algorithms are suitable for VLSI implementation and require only a slight modification of the original CORDIC algorithm.>

Fast evaluation of polynomials is a major goal of computer science, since any continuous function may be approximated as accurately as desired by a polynomial. For instance, most part of current computers evaluate elementary functions using polynomial or rational approximations. J. Duprat and J.M. Muller (1988) presented a new operator, a polynomier, suitable for VLSI implementation, and specifically designed for polynomials computations. This polynomier is composed of two pipe-lined subparts : a squarer (i.e. an operator able to compute the square of a number), and a binomier (i.e. an element which computes expressions of the form Ax+B). The authors present some possible applications of the polynomier. They recall briefly the main characteristics of the architecture. They describe a VLSI implementation of this architecture and present an extension of this architecture to the two's complement computation.>

Complexity results that allow the exact determination or bounding of the online delay of most common arithmetic and elementary functions are presented. These results show that many classical online operators presented in the literature are optimal in delay (but not necessarily in period). The authors propose a way to conserve, for large numbers of manipulations, the main advantage of online arithmetic (the capability of digit-level pipelining) by presenting sparse online arithmetic.>

The authors deal with the detailed VLSI implementation of a fast bit-serial operator designed to perform very high precision (600 decimal digits) additions, multiplications, and divisions, and some of the applications of the circuit are discussed. Online arithmetic needs carry-free redundant number systems. Frequently, the radix chosen is different from 2, since a carry-free addition algorithm can be used in radix r not=2. In radix 2, carry-free addition is possible, but with two inconveniences: the algorithm seems more complicated, and the delay is larger. The authors show that the first inconvenience vanishes if good binary representation of the digits in radix-2 signed digit notation is chosen.>

We describe a general VLSI architecture for the computation of arithmetic expressions including floating-point trancendental functions. This architecture is divided in three parts: a communication machine, the control part of a computation machine and the operative part of this computation machine. In order to compute the most usual trancendental functions, we introduced some general algorithms, presented briefly here, including as a particular case the CORDIC scheme. Our major architecture goals were regularity, parametrization and automatic design. The final chip is designed in a 2-Alu CMOS technology, and its name is FELIN (“Fonctions ELémentaires INtégrées is the french for integrated elementary functions”). This work was supported in part by the GRECO C3and the GCIS of the French CNRS.

In this paper, we present a way to obtain efficient carry-skip adders, built with blocks of different sizes in VLSI technologies. We give some results about two-level carry-skip adders. We reduce our optimization problem to a geometrical problem, solved by means of an algorithm easily implemented on a microcomputer. Then we present an example of the realization of such an adder.

We give necessary and sufficient conditions in order that the infinite product or sum of the terms of a positive decreasing sequence generates the reals in a given interval.