Abstract: RISC-V is a promising free and open-source instruction set architecture. Most of the instruction set has been standardized and several hardware implementations are commercially available. In this paper we highlight features of RISC-V that are interesting for optimizing implementations of cryptographic primitives. We provide the first optimized assembly implementations of table-based AES, bitsliced AES, ChaCha, and the Keccak-\(f\)[1600] permutation for the RV32I instruction set. With respect to public-key cryptography, we study the performance of arbitrary-precision integer arithmetic without a carry flag. We then estimate the improvement that can be gained by several RISC-V extensions. These performance studies also serve to aid design choices for future RISC-V extensions and implementations.

Abstract: The most popular algorithms for computing the matrix exponential are those based on the scaling and squaring technique. For optimal efficiency these are usually tuned to a particular precision of floating-point arithmetic. We design a new scaling and squaring algorithm that takes the unit roundoff of the arithmetic as input and chooses the algorithmic parameters in order to keep the forward error in the underlying Pade approximation below the unit roundoff. To do so, we derive an explicit expression for all the coefficients in an error expansion for Pade approximants to the exponential and use it to obtain a new bound for the truncation error. We also derive a new technique for selecting the internal parameters used by the algorithm, which at each step decides whether to scale or to increase the degree of the approximant. The algorithm can employ diagonal Pade approximants or Taylor approximants and can be used with a Schur decomposition or in transformation-free form. Our numerical experiments show that the new algorithm performs in a forward stable way for a wide range of precisions and that the most accurate of our implementations, the Taylor-based transformation-free variant, is superior to existing alternatives.

Abstract: This paper proposes an innovative Floating Point (FP) architecture for Variable Precision (VP) computation suitable for high precision FP computing, based on a refined version of the UNUM type I format. This architecture supports VP FP intervals where each interval endpoint can have up to 512 bits of mantissa. The proposed hardware architecture is pipelined and has an internal word-size of 64 bits. Computations on longer mantissas are performed iteratively on the existing hardware. The prototype is integrated in a RISC-V environment, it is exposed to the user through an instruction set extension. The paper we provide an example of software usage. The system has been prototyped on a FPGA (Field-Programmable Gate Array) platform and also synthesized for a 28nm FDSOI process technology. The respective working frequency of FPGA and ASIC implementations are 50MHz and 600MHz. The estimated chip area is 1.5mm2 and the estimated power consumption is 95mW. The flops performance of this architecture remains within the range of a regular fixed-precision IEEE FPU while enabling arbitrary precision computation at reasonable cost.

Abstract: We present algorithms for real and complex dot product and matrix multiplication in arbitrary-precision floating-point and ball arithmetic. A low-overhead dot product is implemented on the level of GMP limb arrays; it is about twice as fast as previous code in MPFR and Arb at precision up to several hundred bits. Up to 128 bits, it is 3-4 times as fast, costing 20-30 cycles per term for floating-point evaluation and 40-50 cycles per term for balls. We handle large matrix multiplications even more efficiently via blocks of scaled integer matrices. The new methods are implemented in Arb and significantly speed up polynomial operations and linear algebra.

Abstract: In this study, we implemented large integer multiplication with Single Instruction Multiple Data (SIMD) instructions. We evaluated the implementation on a processor with Cannon Lake microarchitecture, containing Intel AVX-512IFMA (Integer Fused Multiply-Add) instructions. AVX-512IFMA can compute multiple 52-bit integer multiplication and addition operations through one instruction and it has the potential to process large integer multiplications faster than its conventional AVX-512 counterpart. Furthermore, the AVX-512IFMA instructions take three 52-bit integers of 64-bit spaces as operands, and we can use the remaining 12 bits effectively to accumulate carries (reduced-radix representation). For multiplication in the context of larger integers, we applied the Karatsuba and Basecase methods to our program. The former is known to be a faster algorithm than the latter. For evaluation purposes, we compared execution times against extant alternatives and the GNU Multiple Precision Arithmetic Library (GMP). This comparison showed that we were able to achieve a substantive improvement in performance. Specifically, our proposed approach was up to approximately 3.07 times faster than AVX-512F (Foundation) and approximately 2.97 times faster than GMP.

Abstract: Recent studies have shown that High-Level Synthesis (HLS) is an efficient way to design operators for floating-point arithmetic, or for emerging alternative formats such as posits. However, HLS tools support different supersets of different subsets of the C language -- for example, support for arbitrary-sized bit vectors may be provided through vendor-specific data-type libraries such as ac_int, ap_int, or int1 to int64, while others only support the standard C integer types. This is a problem when carefully tuning an operator's internal data-path, as there is no portable HLS standard for arbitrary width integers, and vendor libraries may introduce implicit casts and extensions that can hide subtle bugs. Each vendor also offers varying support for important operator-building primitives, such as platform-optimized leading-zero count. To address such problems, this work introduces Hint (hardware integer), a header-only compatibility layer offering a consistent and comprehensive interface to signed and unsigned arbitrary-sized integers. To avoid bugs Hint is strongly typed, requiring exact matching of expression widths and types -- this type-checking is performed statically using the C++ template system, and adds no overhead at synthesis time. The current implementation wraps ac_int and ap_int with no performance or resource overhead when synthesized on Xilinx or Intel FPGAs. It also offers a Boost::multiprecision backend for fast simulation. Hint is open-source and extensible, and aims to provide an optimized superset of existing library primitives. This work is evaluated with arithmetic operators useful when implementing floating-point and posit operators (shifter, leading zero counter, fused shifter+sticky) deployed using two mainstream HLS tools (Xilinx VivadoHLS, and IntelHLS). A complete posit adder operator has also been written using Hint, showing no overhead when compared to the original operator written for Xilinx FPGAs.

Abstract: We describe algorithms to compute elliptic functions and their relatives (Jacobi theta functions, modular forms, elliptic integrals, and the arithmetic-geometric mean) numerically to arbitrary precision with rigorous error bounds for arbitrary complex variables. Implementations in ball arithmetic are available in the open source Arb library. We discuss the algorithms from a concrete implementation point of view, with focus on performance at tens to thousands of digits of precision.

Abstract: Error-free transformation (EFT) has been recently applied to solve ill-conditioned problems. This transformation can reduce the number of arithmetic operations required compared to multiple precision arithmetic. In this study, we implement double-double (DD) BLAS1 functions with EFT and propose the application of the approach to explicit extrapolation methods for solving initial value problems of ordinary differential equations (ODEs). The presented routines can be effective for a large system of linear ODEs, especially when a harmonic sequence is used.

Abstract: The Vold–Kalman (VK) order tracking filter plays a vital role in the order analysis of noise in various fields. However, owing to the limited accuracy of double-precision floating-point data type, the order of the filter cannot be too high. This problem of accuracy makes it impossible for the filter to use a smaller bandwidth, meaning that the extracted order signal has greater noise. In this paper, the Python mpmath arbitrary-precision floating-point arithmetic library is used to implement a high-order VK filter. Based on this library, a filter with arbitrary bandwidth and arbitrary difference order can be implemented whenever necessary. Using the proposed algorithm, a narrower transition band and a flatter passband can be obtained, a good filtering effect can still be obtained when the sampling rate of the speed signal is far lower than that of the measured signal, and it is possible to extract narrowband signals from signals with large bandwidth. Test cases adopted in this paper show that the proposed algorithm has better filtering effect, better frequency selectivity, and stronger anti-interference ability compared with double-precision data type algorithm.

Abstract: In this paper, we present a digital circuit of arithmetic unit implementing addition and subtraction operations in multiple-precision arithmetic (MPA). This adder-subtractor unit is a part of MPA coprocessor supporting and offloading the central processing unit (CPU) in computations requiring precision higher than 32/64 bits. Although addition and subtraction operations of two n-digit numbers require O(n) operations, the efficient implementation of these operations can provide valuable time-savings for the MPA coprocessor. Furthermore, MPA numbers are usually stored with the use of the sign-magnitude representation which is not so straightforward for addition/subtraction implementation as the two’s complement representation.Our adder-subtractor unit is implemented using the very high speed integrated circuit hardware description language (VHDL) and benchmarked on Xilinx Artix-7 FPGA. The developed digital circuit of the MPA adder-subtractor works with integer numbers of precision varying in the range between 64 bits and 32 kbits with the limb size set to 64 bits. It can currently work with the clock frequency exceeding 450 MHz. For the developed implementation, the addition of two k-limb numbers takes 33+k clock cycles. Hence, the developed coprocessor is 1.7 times faster than a single core of modern i7 processor for precision set to 32704 bits.

Abstract: Many algorithms feature an iterative loop that converges to the result of interest. The numerical operations in such algorithms are generally implemented using finite-precision arithmetic, either fixed- or floating-point, most of which operate least-significant digit first. This results in a fundamental problem: if, after some time, the result has not converged, is this because we have not run the algorithm for enough iterations or because the arithmetic in some iterations was insufficiently precise? There is no easy way to answer this question, so users will often over-budget precision in the hope that the answer will always be to run for a few more iterations. We propose a fundamentally new approach: with the appropriate arithmetic able to generate results from most-significant digit first, we show that fixed compute-area hardware can be used to calculate an arbitrary number of algorithmic iterations to arbitrary precision, with both precision and approximant index increasing in lockstep. Consequently, datapaths constructed following our principles demonstrate efficiency over their traditional arithmetic equivalents where the latter's precisions are either under- or over-budgeted for the computation of a result to a particular accuracy. Use of most-significant digit-first arithmetic additionally allows us to declare certain digits to be stable at runtime, avoiding their recalculation in subsequent iterations and thereby increasing performance and decreasing memory footprints. Versus arbitrary-precision iterative solvers without the optimisations we detail herein, we achieve up-to 16$\times$ performance speedups and 1.9x memory savings for the evaluated benchmarks.

Abstract: The IEEE 1788–2015 has standardized interval arithmetic. However, few libraries for interval arithmetic are compliant with this standard. In the first part of this paper, the main features of the IEEE 1788–2015 standard are detailed, namely the structure into 4 levels, the possibility to accomodate a new mathematical theory of interval arithmetic through the notion of flavor, and the mechanism of decoration for handling exceptions. These features were not present in the libraries developed prior to the elaboration of the standard. MPFI is such a library: it is a C library, based on MPFR, for arbitrary precision interval arithmetic. MPFI is not (yet) compliant with the IEEE 1788–2015 standard for interval arithmetic: the planned modifications are presented. Some considerations about performance and HPC on interval computations based on this standard, or on MPFI, conclude the paper.

Abstract: Functions of matrices arise in numerous applications, and their accurate and efficient evaluation is an important topic in numerical linear algebra. In this thesis, we explore methods to compute them reliably in arbitrary precision arithmetic: on the one hand, we develop some theoretical tools that are necessary to reduce the impact of the working precision on the algorithmic design stage; on the other, we present new numerical algorithms for the evaluation of primary matrix functions and the solution of matrix equations in arbitrary precision environments. Many state-of-the-art algorithms for functions of matrices rely on polynomial or rational approximation, and reduce the computation of f(A) to the evaluation of a polynomial or rational function at the matrix argument A. Most of the algorithms developed in this thesis are no exception, thus we begin our investigation by revisiting the Paterson-Stockmeyer method, an algorithm that minimizes the number of nonscalar multiplications required to evaluate a polynomial of a certain degree. We introduce the notion of optimal degree for an evaluation scheme, and derive formulae for the sequences of optimal degree for the schemes used in practice to evaluate truncated Taylor and diagonal Pade approximants. If the rational function r approximates f, then it is reasonable to expect that a solution to the matrix equation r(X) = A will approximate the functional inverse of f. In general, infinitely many matrices can satisfy this kind of equation, and we propose a classification of the solutions that is of practical interest from a computational standpoint. We develop a precision-oblivious numerical algorithm to compute all the solutions that are of interest in practice, which behaves in a forward stable fashion. After establishing these general techniques, we concentrate on the matrix exponential and its functional inverse, the matrix logarithm. We present a new scaling and squaring approach for computing the matrix exponential in high precision, which combines a new strategy to choose the algorithmic parameters with a bound on the forward error of Pade approximants to the exponential. Then, we develop two algorithms, based on the inverse scaling and squaring method, for evaluating the matrix logarithm in arbitrary precision. The new algorithms rely on a new forward error bound for Pade approximants, which for highly nonnormal matrices can be considerably smaller than the classic bound of Kenney and Laub. Our experimental results show that in double precision arithmetic the new approaches are comparable with the state-of-the-art algorithm for computing the matrix logarithm, and experiments in higher precision support the conclusion that the new algorithms behave in a forward stable way, typically outperforming existing alternatives. Finally, we consider a problem of the form f(A)b, and focus on methods for computing the action of the weighted geometric mean of two large and sparse positive definite matrices on a vector. We present two new approaches based on numerical…

Abstract: In the present study a new multiple-precision arithmetic environment in MATLAB is developed based on exflib, a fast multiple-precision arithmetic in the programming language C, C++ and FORTRAN. We discuss design and implementation of interface between MATLAB and exflib without any modifications to exflib. Although the proposed design incurs the overhead in execution, it is more accurate and faster than Variable Precision Arithmetic (VPA), which is the official multiple-precision arithmetic environment in MATLAB. We also exhibit an efficiency of the proposed environment for reliable computation to overcome the difficulties caused by rounding errors.

Abstract: In this text explicit forms of several higher precision order kernel functions (to be used in he differentiation-by-integration procedure) are given for several derivative orders. Also, a system of linear equations is formulated which allows to construct kernels with an arbitrary precision for an arbitrary derivative order. A computer study is realized and it is shown that numerical differentiation based on higher precision order kernels performs much better (w.r.t. errors) than the same procedure based on the usual Legendre-polynomial kernels. Presented results may have implications for numerical implementations of the differentiation-by-integration method.

Abstract: This chapter deals with fundamental theories on the accuracy of numerical calculation and some cases that seems to be important, somewhat different from previous chapters. We must remember that numerical errors are included in the output data of the computer. In particular, do not overlook the important points you need to know when parallelizing codes. Pursuit of calculation speed is, of course, the central theme of this book, however, it is premised that it produces correct results. This chapter introduces a numerical computation method with guaranteed accuracy in large-scale numerical computations, convergence accuracy problems in parallel computing, and high-precision calculation in HPC.

Abstract: We present algorithms for real and complex dot product and matrix multiplication in arbitrary-precision floating-point and ball arithmetic. A low-overhead dot product is implemented on the level of GMP limb arrays; it is about twice as fast as previous code in MPFR and Arb at precision up to several hundred bits. Up to 128 bits, it is 3-4 times as fast, costing 20-30 cycles per term for floating-point evaluation and 40-50 cycles per term for balls. We handle large matrix multiplications even more efficiently via blocks of scaled integer matrices. The new methods are implemented in Arb and significantly speed up polynomial operations and linear algebra.

Abstract: The purpose of the work, the results of which are presented in the article, was to develop a computerized mathematical model of the ladder filter to study the features of the synthesis of matching chains of any order on its basis. To achieve this goal, all the synthesis tasks were solved, including the choice of the prototype, frequency transposition, calculation of the poles of the transfer function of the prototype, calculation of the poles of the transfer function of the filter, calculation of the reflection coefficient of the filter, the calculation of the input resistance of the filter, the implementation of the ladder chain by the Cauer method and denormalization of the values of the elements of the ladder chain. Modeling of characteristics of ladder filters of any high order with frequency characteristic of Chebyshev on the computer with 32-bit operating system is carried out. It is shown that the use of standard mathematics leads to a significant increase in the error of matching of loads at increasing the order of the chain. To increase the accuracy of calculations it is proposed to use software that allows implementing mathematical operations on variables with any length of the mantissa. There are computer libraries where numbers are presented in the form of string type variables and arithmetic operations are carried out on them according to school rules “ in a column ” . Once the calculations are completed, the strings are converted back to normal numbers. The article considers the application of the variant of BigNumber library for high level JavaScript language. To estimate the accuracy of the calculation it is proposed to apply the property of antimetric ladder structure. From the presented results of the calculation it follows that in order to obtain not less than 15 reliable digits after the decimal point for the parameters of the n -th order filter it is necessary to increase the length of the variable mantissa to the value 2n .

Abstract: The popularity and community-driven development model of RISC-V have opened many areas of investigation to researchers and engineers. To overcome some of the IEEE 754 standard's limitations, one currently emerging avenue for computer architecture and systems research is the area of alternative floating-point computation. The UNUM format, for instance, offers variable precision and much flexibility useful to scientific computing or computational geometry. Programmers usually rely on arbitrary precision libraries such as MPFR (itself depending on GMP). However, there is currently no specialized RISC-V support for these libraries, and little support for variable precision arithmetic across the tool chain in general. We propose a framework to explore the potential of variable precision arithmetic in scientific computing applications on RISC-V processors. This work comprises: (i) a floating-point RISC-V copro-cessor which improve accuracy using the UNUM format; (ii) an ISA extension of the RISC-V ISA for the unit, (iii) a programming model for this extension, and (iv) RISC-V optimized routines for the GMP library. Comparing our solution with MPFR on linear systems solvers, we are able to achieve speedups of up to 18× while keeping computational errors within the same order of magnitude. For 512 bits of precision, speedup between 9x and 16x are observed.

Abstract: Quadratic optimization problems (QPs) are ubiquitous, and solution algorithms have matured to a reliable technology. However, the precision of solutions is usually limited due to the underlying floating-point operations. This may cause inconveniences when solutions are used for rigorous reasoning. We contribute on three levels to overcome this issue. First, we present a novel refinement algorithm to solve QPs to arbitrary precision. It iteratively solves refined QPs, assuming a floating-point QP solver oracle. We prove linear convergence of residuals and primal errors. Second, we provide an efficient implementation, based on SoPlex and qpOASES that is publicly available in source code. Third, we give precise reference solutions for the Maros and Meszaros benchmark library.

Abstract: Polyhedral projection is a main operation of the polyhedron abstract domain.It can be computed via parametric linear programming (PLP), which is more efficient than the classic Fourier-Motzkin elimination method.In prior work, PLP was done in arbitrary precision rational arithmetic.In this paper, we present an approach where most of the computation is performed in floating-point arithmetic, then exact rational results are reconstructed.We also propose a workaround for a difficulty that plagued previous attempts at using PLP for computations on polyhedra: in general the linear programming problems are degenerate, resulting in redundant computations and geometric descriptions.

Abstract: In many scientific applications, it is necessary to compute the sums of floating-point numbers. Summation is a building block for many numerical algorithms, such as dot product, Taylor series, polynomial interpolation and numerical integration. However, the summation of large sets of numbers in finite-precision IEEE 754 arithmetic can be very inaccurate due to the accumulation of rounding errors. There are various ways to diminish rounding errors in the floating-point summation. One of them is the use of multiple-precision arithmetic libraries. Such libraries provide data structures and subroutines for processing numbers whose precision exceeds the IEEE 754 floating-point formats. In this paper, we consider multiple-precision summation on hybrid CPUGPU platforms using MPRES, a new software library for multiple-precision computations on CPUs and CUDA compatible GPUs. Unlike existing multiple-precision libraries based on the binary representation of numbers, MPRES uses a residue number system (RNS). In RNS, the number is represented as a tuple of residues obtained by dividing this number by a given set of moduli, and multiple-precision operations such as addition, subtraction and multiplication are naturally divided into groups of reduced-precision operations on residues, performed in parallel and without carry propagation. We consider the algorithm for the addition of multiple-precision floatingpoint numbers in MPRES, as well as three summation algorithms: (1) recursive summation, (2) pairwise summation, and (3) block-parallel hybrid CPU-GPU summation. Experiments show that the hybrid algorithm allows the full utilization of the GPU’s resources, and therefore demonstrates better performance. On the other hand, the parallel computation of the digits (residues) of multiple-precision significands in RNS reduces computation time

Abstract: These course notes are about computing modular forms and some of their arithmetic properties. Their aim is to explain and prove the modular symbols algorithm in as elementary and as explicit terms as possible, and to enable the devoted student to implement it over any ring (such that a sufficient linear algebra theory is available in the chosen computer algebra system). The chosen approach is based on group cohomology and along the way the needed tools from homological algebra are provided.