A parallel algorithm for accurate polynomial evaluation is proposed for SIMD architectures. This is a parallelized version of the compensated Horner scheme using error-free transformations. The proposed parallel algorithm in this paper is fast and is designed to achieve a result as if computed in twice the working precision and then rounded to the working precision. Numerical results are presented showing the performance of this new parallel algorithm.

Mutation testing is an established approach for checking whether code satisfies a code-independent functional specification, and for evaluating whether a test set is adequate. Current mutation testing approaches, however, do not account for accuracy requirements that appear with numerical specifications implemented in floating- point arithmetic code, but which are a frequent part of safety-critical software. We present Magneto, an instantiation of mutation testing that fully automatically generates a test set from a real-valued specification. The generated tests check numerical code for accuracy, robustness and functional behavior bugs. Our technique is based on formulating test case and oracle generation as a constraint satisfaction problem over interval domains, which soundly bounds errors, but is nonetheless efficient. We evaluate Magneto on a standard floating-point benchmark set and find that it outperforms a random testing baseline for producing useful adequate test sets.

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.

Deep Neural Networks (DNN) represent a performance-hungry application. Floating-Point (FP) and custom floating-point-like arithmetic satisfies this hunger. While there is need for speed, inference in DNNs does not seem to have any need for precision. Many papers experimentally observe that DNNs can successfully run at almost ridiculously low precision. The aim of this paper is two-fold: first, to shed some theoretical light upon why a DNN's FP accuracy stays high for low FP precision. We observe that the loss of relative accuracy in the convolutional steps is recovered by the activation layers, which are extremely well-conditioned. We give an interpretation for the link between precision and accuracy in DNNs. Second, the paper presents a software framework for semi-automatic FP error analysis for the inference phase of deep-learning. Compatible with common Tensorflow/Keras models, it leverages the frugally-deep Python/C++ library to transform a neural network into C++ code in order to analyze the network's need for precision. This rigorous analysis is based an Interval and Affine arithmetics to compute absolute and relative error bounds for a DNN. We demonstrate our tool with several examples.

In this paper we target the Fixed-Point (FxP) implementation of Linear Time-Invariant (LTI) filters evaluated with statespace equations. We assume that wordlengths are fixed and that our goal is to determine binary point positions that guarantee the absence of overflows while maximizing accuracy. We provide a model for the worst-case error analysis of FxP filters that gives tight bounds on the output error. Then we develop an algorithm for the determination of binary point positions that takes rounding errors and their amplification fully into account. The proposed techniques are rigorous, i.e., based on proofs, and no simulations are ever used. In practice, Floating-Point (FP) errors that occur in the implementation of FxP design routines can lead to overestimation/underestimation of resulting parameters. Thus, along with FxP analysis of digital filters, we provide FP analysis of our filter design algorithms. In particular, the core measure in our approach, Worst-Case Peak Gain, is defined as an infinite sum and has matrix powers in it. We provide fine-grained FP error analysis of its evaluation and develop multiple precision algorithms that dynamically adapt their internal precision to satisfy an a priori absolute error bound. Our techniques on multiple precision matrix algorithms, such as eigendecomposition, are of independent interest as a contribution to Computer Arithmetic. All algorithms are implemented as C libraries, integrated into an open-source filter code generator and tested on numerical examples.

Polynomial evaluation is a critical part of the efficient floating-point approximation of elementary functions, in software as well as in FPGA-based systems. Designing an optimized polynomial evaluation scheme is a complex and tedious task, due to multitudes of choices in numerous dimensions: the evaluation scheme, like Horner or Estrin, needs to be selected based on implementation goals (latency, throughput, accuracy. . . ) and be adapted to a given architecture, for example by adapting the level of parallelism to the architecture capabilities. For each operation, a fixed-point or floating-point format needs to be chosen, e.g. between formats such as binary32, binary64. Furthermore some schemes and formats induce compromises, in particular when it comes to vectorized evaluation schemes. As part of a longer automated code generation toolchain, polynomial evaluation gains to be used repeatedly. Several aspects of polynomial evaluation have been presented before, such as code generation for Horner schemes with floating-point expansions or optimization of polynomial evaluation schemes. In this work we study both combination and extension of these techniques, striving for their integration in a code generator. In particular, we present an algorithm within the Metalibm-ludgdunum code generation framework, based on input by Metalibm-lutetia. Our intent is to offer state of the art multi-word evaluation with polynomial scheme space exploration with CGPE, Gappa correctness proof and advanced code generation, suited for High-Level Synthesis.

The IEEE 754-2008 Standard governs Floating-Point Arithmetic in all types of Computer Systems. The Standard provides for two radices, 2 and 10. It specifies conversion operations between these radices, but does not allow floating-point formats of different radices to be mixed in computational operations. In contrast, the Standard does provide for mixing formats of one radix in one operation. In order to enhance the Standard and make it closed under all basic computational operations, we propose an algorithm for a correctly rounded mixed-radix Fused-Multiply-and-Add (FMA). Our algorithm takes any combination of IEEE754 binary64 and decimal64 numbers in argument and provides a result in IEEE754 binary64 and decimal64, rounded according to any for the five IEEE754 rounding modes. Our implementation does not require any dynamic memory allocation; its runtime can be bounded statically. We compare our implementation to a basic mixed-radix FMA implementation based on the GMP Multiple Precision library.

Reliable implementation of digital filters in finiteprecision is based on accurate error analysis. However, a small error in the time domain does not guarantee that the implemented filter verifies the initial band specifications in the frequency domain. We propose a novel certified algorithm for the verification of a filter's transfer function, or of an existing finite-precision implementation. We show that this problem boils down to the verification of bounds on a rational function, and further to the positivity of a polynomial. Our algorithm has reasonable runtime efficiency to be used as a criterion in large implementation space explorations. We ensure that there are no false positives but false negative answers may occur. For negative answers we give a tight bound on the margin of acceptable specifications.We demonstrate application of our algorithm to the comparison of various finite-precision implementations of filters already fully designed.

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.

A typical floating-point environment includes support for a small set of about 30 mathematical functions such as exponential, logarithm, trigonometric and hyperbolic functions. These functions are provided by mathematical software libraries (libm), typically in IEEE754 single, double and quad precision. This article suggests to replace this libm paradigm by a more general approach: the on-demand generation of numerical function code, on arbitrary domains and with arbitrary accuracies. First, such code generation opens up the libm function space available to programmers. It may capture a much wider set of functions, and may capture even standard functions on non-standard domains and accuracy/performance points. Second, writing libm code requires fine-tuned instruction selection and scheduling for performance, and sophisticated floating-point techniques for accuracy. Automating this task through code generation improves confidence in the code while enabling better design space exploration, and therefore better time to market, even for the libm functions. This article discusses the new challenges of this paradigm shift, and presents the current state of open-source function code generators available on http://www.metalibm.org/.

This work introduces an approach to the computer-assisted implementation of mathematical functions geared toward special functions such as those occurring in mathematical physics. The general idea is to start with an exact symbolic representation of a function and automate as much as possible of the process of implementing it. In order to deal with a large class of special functions, our symbolic representation is an implicit one: the input is a linear differential equation with polynomial coefficients along with initial values. The output is a C program to evaluate the solution of the equation using domain splitting, argument reduction and polynomial approximations in double-precision arithmetic, in the usual style of mathematical libraries. Our generation method combines symbolic-numeric manipulations of linear ODEs with interval-based tools for the floating-point implementation of "black-box" functions. We describe a prototype code generator that can automatically produce implementations on moderately large intervals. Implementations on the whole real line are possible in some cases but require manual tool setup and code integration. Due to this limitation and as some heuristics remain, we refer to our method as "semi-automatic" at this stage. Along with other examples, we present an implementation of the Voigt profile with fixed parameters that may be of independent interest.

The worst-case peak gain (WCPG) of a linear filter is an important measure for the implementation of signal processing algorithms. It is used in the error propagation analysis for filters, thus a reliable evaluation with controlled precision is required. The WCPG is computed as an infinite sum and has matrix powers in each summand. We propose a direct formula for the lower bound on truncation order of the infinite sum in dependency of desired truncation error. Several multiprecision methods for complex matrix operations are developed and their error analysis performed. A multiprecision matrix powering method is presented. All methods yield a rigorous solution with an absolute error bounded by an a priori given value. The results are illustrated with numerical examples.

In this article, we present an efficient algorithm to compute the faithful rounding of the l 2 -norm of a floating-point vector. This means that the result is accurate to within 1 bit of the underlying floating-point type. This algorithm does not generate overflows or underflows spuriously, but does so when the final result calls for such a numerical exception to be raised. Moreover, the algorithm is well suited for parallel implementation and vectorization. The implementation runs up to 3 times faster than the netlib version on current processors.

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.

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

High confidence in floating-point programs requires proving numerical properties of final and intermediate values. One may need to guarantee that a value stays within some range, or that the error relative to some ideal value is well bounded. This certification may require a time-consuming proof for each line of code, and it is usually broken by the smallest change to the code, e.g., for maintenance or optimization purpose. Certifying floating-point programs by hand is, therefore, very tedious and error-prone. The Gappa proof assistant is designed to make this task both easier and more secure, due to the following novel features: It automates the evaluation and propagation of rounding errors using interval arithmetic. Its input format is very close to the actual code to validate. It can be used incrementally to prove complex mathematical properties pertaining to the code. It generates a formal proof of the results, which can be checked independently by a lower level proof assistant like Coq. Yet it does not require any specific knowledge about automatic theorem proving, and thus, is accessible to a wide community. This paper demonstrates the practical use of this tool for a widely used class of floating-point programs: implementations of elementary functions in a mathematical library.

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.

In many numerical programs there is a need for a high-quality floating-point approximation of useful functions f, such as such as exp, sin, erf. In the actual implementation, the function is replaced by a polynomial p, which leads to an approximation error (absolute or relative) epsiv = p - s or epsiv = p/f -1. The tight yet certain bounding of this error is an important step towards safe implementations.The problem is difficult mainly because that approximation error is very small and the difference p-f is subject to high cancellation. Previous approaches for computing the supremum norm in this degenerate case, have proven to be unsafe, not sufficiently tight or too tedious in manual work.We present a safe and fast algorithm that computes a tight lower and upper bound for the supremum norms of approximation errors. The algorithm is based on a combination of several techniques, including enhanced interval arithmetic, automatic differentiation and isolation of the roots of a polynomial. We have implemented our algorithm and give timings on several examples.

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.