Floating-point Approximation Research Articles

Floating-Point Approximation Enabling Cost-Effective and High-Precision Digital Implementation of FitzHugh-Nagumo Neural Networks.

The study of neuron interactions and hardware implementations are crucial research directions in neuroscience, particularly in developing large-scale biological neural networks. The FitzHugh-Nagumo (FHN) model is a popular neuron model with highly biological plausibility, but its complexity makes it difficult to apply at scale. This paper presents a cost-saving and improved precision approximation algorithm for the digital implementation of the FHN model. By converting the computational data into floating-point numbers, the original multiplication calculations are replaced by adding the floating-point exponent part and fitting the mantissa part with piecewise linear. In the hardware implementation, shifters and adders are used, greatly reducing resource overhead. Implementing FHN neurons by this approximation calculations on FPGA reduces the normalized root mean square error (RMSE) to 3.5% of the state-of-the-art (SOTA) while maintaining a performance overhead ratio improvement of 1.09 times. Compared to implementations based on approximate multipliers, the proposed method achieves a 20% reduction in error at the cost of a 2.8% increase in overhead.This model gained additional biological properties compared to LIF while reducing the deployment scale by only 9%. Furthermore, the hardware implementation of nine coupled circular networks with eight nodes and directional diffusion was carried out to demonstrate the algorithm's effectiveness on neural networks. The error decreased to 60% compared to the single neuron of the SOTA. This hardware-friendly algorithm allows for the low-cost implementation of high-precision hardware simulation, providing a novel perspective for studying large-scale, biologically plausible neural networks.

Algebraic number fields and the LLL algorithm

In this paper we analyze the computational costs of various operations and algorithms in algebraic number fields using exact arithmetic. Let K be an algebraic number field. In the first half of the paper, we calculate the running time and the size of the output of many operations in K in terms of the size of the input and the parameters of K. We include some earlier results about these, but we go further than them, e.g. we also analyze some R-specific operations in K like less-than comparison.In the second half of the paper, we analyze two algorithms: the Bareiss algorithm, which is an integer-preserving version of the Gaussian elimination, and the LLL algorithm, which is for lattice basis reduction. In both cases, we extend the algorithm from Zn to Kn, and give a polynomial upper bound on the running time when the computations in K are performed exactly (as opposed to floating-point approximations).

Using the Interval-Symbol Method with Zero Rewriting to Factor Polynomials over Algebraic Number Fields

The ISZ method (Interval-Symbol method with Zero rewriting) based on stabilization theory was proposed to reduce the amount of exact computations as much as possible but obtain the exact results by aid of floating-point computations. In this paper, we applied the ISZ method to Trager's algorithm which factors univariate polynomials over algebraic number fields. By Maple experiments, we show the efficiency of the ISZ method over the purely exact approach which uses exact computations throughout the execution of the algorithm. Furthermore, we propose a new method called the ISZ* method, which is similar to the ISZ method but beforehand excludes insufficient precisions of floating-point approximation by checking the correctness of the obtained supports. We confirmed that the ISZ* method is more effective than the ISZ method when the initially set precision is not sufficiently high.

Toward Accurate and Fast Summation

We introduce a new accurate summation algorithm based on the error-free summation into floating-point buckets. Our algorithm exploits ideas from Zhu and Hayes’ OnlineExactSum , but it uses a significantly smaller number of accumulators and has a better instruction-level parallelism. In the default setting, our implementation aaaSum returns a faithfully rounded floating-point approximation of the true sum. We also discuss possible modifications for the computation of reproducible, correctly rounded, and multiple precision floating-point approximations. The computational overhead for any of these modifications is kept comparably small. Numerical tests demonstrate that aaaSum performs well for very small to large problem sizes, independent of the condition number of the problem. We compare our algorithm with other accurate and high-precision summation approaches.

Indirect Predicates for Geometric Constructions

Geometric predicates are a basic ingredient to implement a vast range of algorithms in computational geometry. Modern implementations employ floating point filtering techniques to combine efficiency and robustness, and state-of-the-art predicates are guaranteed to be always exact while being only slightly slower than corresponding (inexact) floating point implementations. Unfortunately, if the input to these predicates is an intermediate construction of an algorithm, its floating point representation may be affected by an approximation error, and correctness is no longer guaranteed.This paper introduces the concept of indirect geometric predicate: instead of taking the intermediate construction as an explicit input, an indirect predicate considers the primitive geometric elements which are combined to produce such a construction. This makes it possible to keep track of the floating point approximation, and thus to exploit efficient filters and expansion arithmetic to exactly resolve the predicate with minimal overhead with respect to a naive floating point implementation. As a representative example, we show how to extend standard predicates to the case of points of intersection of linear elements (i.e. lines and planes) and show that, on classical problems, this approach outperforms state-of-the-art solutions based on lazy exact intermediate representations.

A note on Dekker\u2019s FastTwoSum algorithm

More than 45 years ago, Dekker proved that it is possible to evaluate the exact error of a floating-point sum with only two additional floating-point operations, provided certain conditions are met. Today the respective algorithm for transforming a sum into its floating-point approximation and the corresponding error is widely referred to as {{,mathrm{FastTwoSum},}}. Besides some assumptions on the floating-point system itself—all of which are satisfied by any binary IEEE 754 standard conform arithmetic, the main practical limitation of {{,mathrm{FastTwoSum},}} is that the summands have to be ordered according to their exponents. In most preceding applications of {{,mathrm{FastTwoSum},}}, however, a more stringent condition is used, namely that the summands have to be sorted according to their absolute value. In remembrance of Dekker’s work, this note reminds the original assumptions for an error-free transformation via{{,mathrm{FastTwoSum},}}. Moreover, we generalize the conditions for arbitrary bases and discuss a possible modification of the {{,mathrm{FastTwoSum},}} algorithm to extend its applicability even further. Subsequently, a range of programs exploiting the wider applicability is presented. This comprises the OnlineExactSum algorithm by Zhu and Hayes, an error-free transformation from a product of three floating-point numbers to a sum of the same number of addends, and an algorithm for accurate summation proposed by Demmel and Hida.

Computing Hypergeometric Functions Rigorously

We present an efficient implementation of hypergeometric functions in arbitrary-precision interval arithmetic. The functions 0 F 1 , 1 F 1 , 2 F 1 , and 2 F 0 (or the Kummer U -function) are supported for unrestricted complex parameters and argument, and, by extension, we cover exponential and trigonometric integrals, error functions, Fresnel integrals, incomplete gamma and beta functions, Bessel functions, Airy functions, Legendre functions, Jacobi polynomials, complete elliptic integrals, and other special functions. The output can be used directly for interval computations or to generate provably correct floating-point approximations in any format. Performance is competitive with earlier arbitrary-precision software and sometimes orders of magnitude faster. We also partially cover the generalized hypergeometric function p F q and computation of high-order parameter derivatives.

Stack and register complexity of radix conversions

We investigate the question of computational resources (such as stacks and counters) necessary to perform radix conversions. To this end it is shown that no PDA can compute the significand of the best n-digit floating point approximation of a power of an incommensurable radix. This extends the results of W. Clinger. We also prove that a two counter machine with input is capable of such conversions. On the other hand we note a curious asymmetry with respect to the order in which the digits are input by showing that a two counter machine can decode its input online if the digits are presented in the most-to-least significant order while no such machine can decode its input in this manner if the digits are presented in the least-to-most significant order. Some structural results about two counter machines (with input) are also established.

Computing Class Polynomials for Abelian Surfaces

We describe a quasilinear algorithm for computing Igusa class polynomials of Jacobians of genus-2 curves via complex floating-point approximations of their roots. After providing an explicit treatment of the computations in quartic CM fields and their Galois closures, we pursue an approach due to Dupont for evaluating ϑ-constants in quasilinear time using Newton iterations on the Borchardt mean. We report on experiments with our implementation and present an example with class number 20 016.

Improved Error Bounds for Inner Products in Floating-Point Arithmetic

Given two floating-point vectors $x,y$ of dimension $n$ and assuming rounding to nearest, we show that if no underflow or overflow occurs, any evaluation order for an inner product returns a floating-point number ${\widehat r}$ such that $|{\widehat r}- x^Ty| \leqslant nu|x|^T|y|$ with $u$ the unit roundoff. This result, which holds for any radix and with no restriction on $n$, can be seen as a generalization of a similar bound given in [S. M. Rump, BIT, 52 (2012), pp. 201--220] for recursive summation in radix $2$, namely, $|{\widehat r}- x^Te| \leqslant (n-1)u|x|^Te$ with $e=[1,1,\ldots,1]^T$. As a direct consequence, the error bound for the floating-point approximation $\widehat{C}$ of classical matrix multiplication with inner dimension $n$ simplifies to $|\widehat{C}-AB|\leqslant nu|A||B|$.

Perturbation Analysis of the QR factor R in the context of LLL lattice basis reduction

In 1982, Arjen Lenstra, Hendrik Lenstra Jr. and László Lovász introduced an efficiently computable notion of reduction of basis of a Euclidean lattice that is now commonly referred to as LLL-reduction. The precise definition involves the R-factor of the QR factorization of the basis matrix. In order to circumvent the use of rational/exact arithmetic with large bit-sizes, it is tempting to consider using floating-point arithmetic with small precision to compute the R-factor. In the present article, we investigate the accuracy of the factor R of the QR factorisation of an LLL-reduced basis. Our main contribution is the first fully rigorous perturbation analysis of the R-factor of LLL-reduced matrices under columnwise perturbations. Our results are very useful to devise LLL-type algorithms relying on floating-point approximations.

Fundamental domains for genus-zero and genus-one congruence subgroups

AbstractIn this paper, we compute Ford fundamental domains for all genus-zero and genus-one congruence subgroups. This is a continuation of previous work, which found all such groups, including ones that are not subgroups ofPSL(2,ℤ). To compute these fundamental domains, an algorithm is given that takes the following as its input: a positive square-free integerf, which determines a maximal discrete subgroup Γ0(f)+ofSL(2,ℝ); a decision procedure to determine whether a given element of Γ0(f)+is in a subgroupG; and the index ofGin Γ0(f)+. The output consists of: a fundamental domain forG, a finite set of bounding isometric circles; the cycles of the vertices of this fundamental domain; and a set of generators ofG. The algorithm avoids the use of floating-point approximations. It applies, in principle, to any group commensurable with the modular group. Included as appendices are: MAGMA source code implementing the algorithm; data files, computed in a previous paper, which are used as input to compute the fundamental domains; the data computed by the algorithm for each of the congruence subgroups of genus zero and genus one; and an example, which computes the fundamental domain of a non-congruence subgroup.

The complexity of class polynomial computation via floating point approximations

We analyse the complexity of computing class polynomials, that are an important ingredient for CM constructions of elliptic curves, via complex floating point approximations of their roots. The heart of the algorithm is the evaluation of modular functions in several arguments. The fastest one of the presented approaches uses a technique devised by Dupont to evaluate modular functions by Newton iterations on an expression involving the arithmetic-geometric mean. It runs in time $O (|D| \log^5 |D| \log \log |D|) = O (|D|^{1 + \epsilon}) = O ( h^{2 + \epsilon})$ for any $\epsilon > 0$, where $D$ is the CM discriminant and $h$ is the degree of the class polynomial. Another fast algorithm uses multipoint evaluation techniques known from symbolic computation; its asymptotic complexity is worse by a factor of $\log |D|$. Up to logarithmic factors, this running time matches the size of the constructed polynomials. The estimate also relies on a new result concerning the complexity of enumerating the class group of an imaginary-quadratic order and on a rigorously proven upper bound for the height of class polynomials.

How to read floating point numbers accurately

Converting decimal scientific notation into binary floating point is nontrivial, but this conversion can be performed with the best possible accuracy without sacrificing efficiency.Consider the problem of converting decimal scientific notation for a number into the best binary floating point approximation to that number, for some fixed precision. This problem cannot be solved using arithmetic of any fixed precision. Hence the IEEE Standard for Binary Floating-Point Arithmetic does not require the result of such a conversion to be the best approximation.This paper presents an efficient algorithm that always finds the best approximation. The algorithm uses a few extra bits of precision to compute an IEEE-conforming approximation while testing an intermediate result to determine whether the approximation could be other than the best. If the approximation might not be the best, then the best approximation is determined by a few simple operations on multiple-precision integers, where the precision is determined by the input. When using 64 bits of precision to compute IEEE double precision results, the algorithm avoids higher-precision arithmetic over 99% of the time.The input problem considered by this paper is the inverse of an output problem considered by Steele and White: Given a binary floating point number, print a correctly rounded decimal representation of it using the smallest number of digits that will allow the number to be read without loss of accuracy. The Steele and White algorithm assumes that the input problem is solved; an imperfect solution to the input problem, as allowed by the IEEE standard and ubiquitous in current practice, defeats the purpose of their algorithm.

Parametrized Family of 2-D Non-factorable FIR Lossless Systems and Gröbner Bases

The factorability of one-dimensional (1-D) FIR lossless transfer matrices [1] in terms of Givens rotations produces the parameters that can be used for an optimal design of filter banks with prespecified filtering characteristics. Two dimensional (2-D) FIR lossless systems behave quite differently, however. Venkataraman-Levy [2] and Basu-Choi-Chiang [3] have constructed 2-D FIR paraunitary matrices of McMillan degrees (2,2) that are not factorable. Because of the state-space realization used in the construction, they are floating-point approximations, and they do not produce explicit parametrizations that can be used for optimal design process. In this paper, we formulate the lossless condition and nonfactorability condition of a 2-D FIR paraunitary matrix using multivariate polynomials in the coefficients. The resulting polynomial system can be explicitly solved with Grobner bases. By studying the polynomial system, we obtain a continuous one parameter family of 2-D 2×2 non-factorable paraunitary matrices. As an example, we get a closed-form expression for a 2-D 2×2 paraunitary matrix that is not factorable into rotations and delays.

Shortest Path Geometric Rounding

Exact implementations of algorithms of computational geometry are subject to exponential growth in running time and space. In particular, coordinate bit-complexity can grow exponentially when algorithms are cascaded : the output of one algorithm becomes the input to the next. Cascading is a significant problem in practice. We propose a geometric rounding technique: shortest path rounding . Shortest path rounding trades accuracy for space and time and eliminates the exponential cost introduced by cascading. It can be applied to all algorithms which operate on planar polygonal regions, for example, set operations, transformations, convex hull, triangulation, and Minkowski sum. Unlike other geometric rounding techniques, shortest path rounding can round vertices to arbitrary lattices, even in polar coordinates, as long as the rounding cells are connected. (Other rounding techniques can only round to the integer grid.) On the integer grid, shortest path rounding introduces less combinatorial change and geometric error than the other rounding methods. Three algorithms are given for shortest path rounding, one of which we have used in industrial application software since 1992. In combination with recent advances in exact floating point evaluation of numerical primitives, shortest path geometric rounding yields a practical solution to numerical issues in computational geometry. Geometric algorithms can be implemented exactly on floating point input coordinates; the exact output coordinates can be rounded to accurate floating point approximations; and the cost of each arithmetic operation is only a little more than if it were implemented as a single hardware floating point operation.

Type matching, type-graphs, and the Schanuel conjecture

This work considers type systems that are defined by type-graphs (tgraphs), which are rooted directed graphs with order among the edges leaving each node. Tgraphs are uniquely mapped into polynomials which, in turn, are each evaluated at a special point to yield an irrational number named the tgraph's magic number . This special point is chosen using the Schanuel conjecture. It is shown that each tgraph can be uniquely represented by this magic number; namely, types are equal if and only if the corresponding magic numbers are equal. Since irrational numbers require infinite precision, the algorithm for generating magic numbers is carried out using a double-precision floating-point approximation. This approximation is viewed as a hashing scheme, mapping the infinite domain of the irrational numbers into finite computer words. The proposed hashing scheme was investigated experimentally, with the conclusion that it is a good and practical hashing method. In tests involving over a million randomly chosen tgraphs, we have not encountered a single collision. We conclude that this method for representation and management of types is practical, and offers novel possibilities for enforcing strict type matching at link time among separately compiled modules.

How to read floating point numbers accurately

Consider the problem of converting decimal scientific notation for a number into the best binary floating point approximation to that number, for some fixed precision. This problem cannot be solved using arithmetic of any fixed precision. Hence the IEEE Standard for Binary Floating-Point Arithmetic does not require the result of such a conversion to be the best approximation. This paper presents an efficient algorithm that always finds the best approximation. The algorithm uses a few extra bits of precision to compute an IEEE-conforming approximation while testing an intermediate result to determine whether the approximation could be other than the best. If the approximation might not be the best, then the best approximation is determined by a few simple operations on multiple-precision integers, where the precision is determined by the input. When using 64 bits of precision to compute IEEE double precision results, the algorithm avoids higher-precision arithmetic over 99% of the time. The input problem considered by this paper is the inverse of an output problem considered by Steele and White: Given a binary floating point number, print a correctly rounded decimal representation of it using the smallest number of digits that will allow the number to be read without loss of accuracy. The Steele and White algorithm assumes that the input problem is solved; an imperfect solution to the input problem, as allowed by the IEEE standard and ubiquitous in current practice, defeats the purpose of their algorithm.

Error-free computer solution of certain systems of linear equations

In this article we propose a procedure which generates the exact solution for the system A x = b , where A is an integral nonsingular matrix and b is an integral vector, by improving the initial floating-point approximation to the solution. This procedure, based on an easily programmed method proposed by Aberth [1], first computes the approximate floating-point solution x * by using an available linear equation solving algorithm. Then it extracts the exact solution x from x * if the error in the approximation x * is sufficiently small. An a posteriori upper bound for the error of x * is derived when Gaussian Elimination with partial pivoting is used. Also, a computable upper bound for |det( A )|, which is an alternative to using Hadamard's inequality, is obtained as a byproduct of the Gaussian Elimination process.

Least significant bit evaluation of arithmetic expressions in single-precision

Single-precision floatingpoint computations may yield an arbitrary false result due to cancellation and rounding errors. This is true even for very simple, structured arithmetic expressions such as Horner's scheme for polynomial evaluation. A simple procedure will be presented for fast calculation of the value of an arithmetic expression to least significant bit accuracy in single precision computation. For this purpose in addition to the floating-point arithmetic only a precise scalar product (cf. [2]) is required. If the initial floatingpoint approximation is not too bad, the computing time of the new algorithm is approximately the same as for usual floating-point computation. If not, the essential progress of the presented algorithm is that the inaccurate approximation is recognized and corrected. The algorithm achieves high accuracy, i.e. between the left and the right bound of the result there is at most one more floating-point number. A rigorous estimation of all rounding errors introduced by floating-point arithmetic is given for general triangular linear systems. The theorem is applied to the evaluation of arithmetic expressions.