We study a bug that we found in the GNU Emacs text editor when built against the Cairo graphics library. We analyze both the Emacs code and the Cairo code, and we suggest what can be done to avoid unexpected results. This involves a particular case with a computation that can be reduced to the equivalent floating-point expression ((1/s) · b) · s, where s and b are small positive integers such that b < s and the basic operations are rounded to nearest. The analysis takes into account the values of s and b that can occur in practice, and the suggestions for workarounds must avoid handling this particular case in a separate branch or breaking the structure of the Cairo library (so that just returning b is not possible).

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

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.

We describe algorithms used to optimize the GNU MPFR library when the operands fit into one or two words. On modern processors, this gives a speedup for a correctly rounded addition, subtraction, multiplication, division or square root in the standard binary64 format (resp. binary128) between 1.8 and 3.5 (resp. between 1.6 and 3.2). We also introduce a new faithful rounding mode, which enables even faster computations. Those optimizations will be available in version 4 of MPFR.

We present a fast algorithm together with its low-level implementation of correctly rounded arbitrary-precision floating-point summation. The arithmetic is the one used by the GNU MPFR library: radix 2; no subnormals; each variable (each input and the output) has its own precision. We also give a worst-case complexity of this algorithm and describe how the implementation is tested.

We present a fast algorithm together with its low-level implementation of correctly rounded arbitrary-precision floating-point summation. The arithmetic is the one used by the GNU MPFR library: radix 2, no subnormals, each variable (each input and the output) has its own precision. We also describe how the implementation is tested.

SIPE (Small Integer Plus Exponent) is a mini-library in the form of a C header file, to perform floating-point computations in very low precisions with correct rounding to nearest in radix 2. The goal of such a tool is to do proofs of algorithms/properties or computations of tight error bounds in these precisions by exhaustive tests, in order to try to generalize them to higher precisions. The currently supported operations are addition, subtraction, multiplication (possibly with the error term), FMA, and miscellaneous comparisons and conversions. Timing comparisons have been done with hardware IEEE-754 floating point and with GNU MPFR.

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.

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.

The correct rounding of the function pow: (x, y) rarrxyis currently based on Ziv's iterative approximation process. In order to ensure its termination, cases when xyfalls on a rounding-boundary must be filtered out. Such rounding-boundaries are floating-point numbers and midpoints between two consecutive floating-point numbers. Detecting rounding-boundaries for pow is a difficult problem. Previous approaches use repeated square root extraction followed by repeated square and multiply. This paper presents a new rounding-boundary test for pow in double precision, which reduces this to a few comparisons with precomputed constants. These constants are deduced from worst cases for the Table Maker's Dilemma, searched over a small subset of the input domain. This is a novel use of such worst-case bounds. The resulting algorithm has been designed for a fast-on-average correctly rounded implementation of pow, considering the scarcity of rounding-boundary cases. It does not stall average computations for rounding-boundary detection. This paper includes its correctness proof and experimental results.

One considers the problem of finding hard to round cases of a periodic function for large floating-point inputs, more precisely when the function cannot be efficiently approximated by a polynomial. This is one of the last few issues that prevents from guaranteeing an efficient computation of correctly rounded transcendentals for the whole IEEE-754 double precision format. The first non-naive algorithm for that problem is presented, with a heuristic complexity of O(20.676p) for a precision of p bits. The efficiency of the algorithm is shown on the largest IEEE-754 double precision binade for the sine function, and some corresponding bad cases are given. We can hope that all the worst cases of the trigonometric functions in their whole domain will be found within a few years, a task that was considered out of reach until now.

This article presents a multiple-precision binary floating-point library, written in the ISO C language, and based on the GNU MP library. Its particularity is to extend to arbitrary-precision, ideas from the IEEE 754 standard, by providing correct rounding and exceptions . We demonstrate how these strong semantics are achieved---with no significant slowdown with respect to other arbitrary-precision tools---and discuss a few applications where such a library can be useful.

This paper presents extensions to Lefevre's algorithm that computes a lower bound on the distance between a segment and a regular grid Zopf2. This algorithm and, in particular, the extensions are useful in the search for worst cases for the exact rounding of unary elementary functions or base-conversion functions. The proof that is presented is simpler and less technical than the original proof. This paper also gives benchmark results with various optimization parameters, explanations of these results, and an application to base conversion

We propose a new algorithm to find worst cases for the correct rounding of a mathematical function of one variable. We first reduce this problem to the real small value problem - i.e., for polynomials with real coefficients. Then, we show that this second problem can be solved efficiently by extending Coppersmith's work on the integer small value problem - for polynomials with integer coefficients - using lattice reduction. For floating-point numbers with a mantissa less than N and a polynomial approximation of degree d, our algorithm finds all worst cases at distance less than N/sup -d2//2d+1 from a machine number in time O(N/sup (d+1/2d+1)+/spl epsiv//). For d=2, a detailed study improves on the O(N/sup 2/(3+/spl epsiv/)/) complexity from Lefevre's algorithm to O(N/sup 4/(7+/spl epsiv/)/). For larger d, our algorithm can be used to check that there exist no worst cases at distance less than N/sup -k/ in time O(N/sup 1/(2+/spl epsiv/)/).

We propose a new algorithm to find worst cases for correct rounding of an analytic function. We first reduce this problem to the real small value problem - i.e. for polynomials with real coefficients. Then we show that this second problem can be solved efficiently, by extending Coppersmith's work on the integer small value problem - for polynomials with integer coefficients - using lattice reduction (D. Coppersmith, 1996; 2001). For floating-point numbers with a mantissa less than N, and a polynomial approximation of degree d, our algorithm finds all worst cases at distance < N/sup -d2//(2d+1) from a machine number in time O(N/sup ((d+1)/(2d+1))+/spl epsiv//). For d=2, this improves on the O(N/sup 2/(3+/spl epsiv/)/) complexity from Lefevre's algorithm (V. Lefevre, 2000; V. Lefevre et al., 2001) to O(N/sup 3/(5+/spl epsiv/)/). We exhibit some new worst cases found using our algorithm, for double-extended and quadruple precision. For larger d, our algorithm can be used to check that there exist no worst cases at distance < N/sup -k/ in time O(N/sup (1/2)+O(1/k)/).

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.

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.

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.