Matei Is

Bogdan Pasca

Fixed-Point Implementations of the Reciprocal, Square Root and Reciprocal Square Root Functions

Implementations of the reciprocal, square root and reciprocal square root often share a common structure. This article is a survey and comparison of methods (with only slight variations for the three cases) for computing these functions. The comparisons are made in the context of the same accuracy target (faithful rounding) and of an arbitrary fixed-point format for the inputs and outputs (precisions of up to 32 bits). Some of the methods discussed might require some form of range reduction, depending on the input's range. The objective of the article is to optimize the use of fixed-size FPGA resources (block multipliers and block RAMs). The discussions and conclusions are based on synthesis results for FPGAs. They try to suggest the best method to compute the previously mentioned fixedpoint functions on a FPGA, given the input precision. This work compares classical methods (direct tabulation, multipartite tables, piecewise polynomials, Taylor-based polynomials, Newton-Raphson iterations). It also studies methods that are novel in this context: the Halley method and, more generally, the Householder method.

I. INTRODUCTION

FPGAs are becoming increasingly large and complex, with recent devices embedding computational features such as single-precision floating-point operators [START_REF] Langhammer | Floating-Point DSP Block Architecture for FPGAs[END_REF]. If two decades ago researchers were customizing floating-point arithmetic in order for single operators to fit entire devices [START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machines[END_REF], the thousands of embedded floating-point adders and multipliers available in the Altera Arria10 devices make possible the rapid construction of highly optimized datapaths which run at FPGA nominal frequencies. Floating-point datapaths can automatically cope with high dynamic data ranges, shortening design time often at the expense of performance. Fixed-point versions of addition and multiplication have traditionally been implemented in FPGAs using dedicated circuitry. Embedded ripple carry adder cells and DSP blocks containing small embedded multipliers allow fixed-point datapaths to provide small implementation footprints and shorter latency. The pitfall when choosing fixed-point arithmetic is that the user needs to manage data ranges manually -a task which requires good algorithm understanding -often regarded as a rare skill.

In a fixed-point datapath, well trimmed adders and multipliers will have efficient implementations. On the other hand, fixed-point function computations (such as √ x, 1/ √ x, 1/x) are often overlooked in terms of resource requirements. Datapaths embedding a mix of these resources need efficient fixedpoint function implementations. This is particularly important since the properties of floating-point arithmetic allows for very efficient implementations of these functions, which may overweigh the advantages of fixed-point adders and multipliers.

This paper focuses on the fixed-point implementation of these functions, for bit-widths ranging from few bits up to 32 bits. The implementations are flexible in terms of input and output formats allowing for fine-grain operator customization. Input range can be restricted to the interval [START_REF] Langhammer | Floating-Point DSP Block Architecture for FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machines[END_REF], or may span the whole input range allowed by the input format. Multiple implementation techniques from the literature are revisited and employed in this generic context, from tabulation, multipartite methods, Taylor series, Newton-Raphson, and higher order methods. A generator which selects the suitable architectures for user-specified constraints was constructed for this paper.

Due to the double-blind review process we are not describing the full details of this generator in this submission.

II. FPGA FEATURES

In this section we briefly present FPGA architectural features relevant to this work. We restrict our study to using Altera FPGA architectures, although the generic nature of the core generator of this paper is actually vendor agnostic.

Out of the 3 FPGA families from Altera we focus on the low-cost CycloneV [START_REF]Cyclone V Device Handbook[END_REF] devices. In a nutshell, FPGA devices are composed of millions of small tables (4-6 inputs) which can be programmed to perform any logic function (of 4-6 inputs). These small tables are connected using a flexible, programmable interconnect network which allows interconnecting any two nodes in the FPGA. In Altera devices, these table based resources are encapsulated in Adaptive Logic Modules (ALMs). Each ALM can be configured as 2 independent 4-input lookup tables (LUT), one 6-input LUT, and many combinations of functions which share inputs. Logic resources are usually reported in terms of ALMs in Altera devices.

Contemporary FPGA devices also contains hundreds (often thousands -depending on the device) small flexible multipliers (packaged in DSP blocks) and memory blocks. The DSP blocks in the CycloneV device can be configured in the modes presented in Table I.

Dedicated memory blocks in CycloneV devices have a capacity of 10Kb each (M10K). The relevant configurations for this works are: 2048x5bits, 1024x10bits, 512x20bits, 256x40bits.

III. BACKGROUND A. Choosing the Right Method

In order to chose the best method for implementing the functions studied in this paper, the criteria on which to judge the different methods first have to be decided. Having multiple criteria has the disadvantage of complicating the selection process, but it also has the advantage of allowing the user to make the right choice on each separate scenario.

The methods presented in the following are differentiated by three main aspects. This short list of criteria is in no way exclusive, but it highlights some of the common design challenges that users are faced with. It also highlights the constant quest for a better use of all the available hardware resources on FPGAs.

Bitwidth: computing just right in terms of precision, means performing meaningful computations and ensuring the target output precision, as well as using a minimal amount of resources. Computing for a specific FPGA target means getting the most out of what the FPGA has to offer. Choosing an implementation method based on the targeted bitwidth can ensure satisfying both points previously made.

Input range: making no assumptions about the inputs can make a design as general as possible. It can also impose additional constraints on the method used (and even make it unsuitable) and incur the use of additional resources. Choosing an implementation while having knowledge about the input ranges allows more specialized solutions.

Approach: table-based, polynomial approximation and iterative methods are the most common classes of methods used for implementing the functions surveyed in this paper. We have provided hybrids between some of these classes. The choice between these approaches is influenced by the previous two criteria, as well as by cost (resource/speed-wise) and the flexibility each method offers.

Before moving on to classifying the different methods into classes, based on the previously mentioned criteria, one important point has to be made. While this article is interested in fixed-point (FxP) implementations, most of the relevant literature deals with floating-point (FP) implementations. This means that previous findings (relevant in FP), may no longer be valid, due to the limited convergence domains of the methods, the need for range reduction, or due to the specificities of the different target platform.

Another significant difference (from most of the existing works) which may bias the findings in this article is allowing the user to specify the input fixed-point format. This is another aspect that is less common in the literature and which justifies to a certain measure the choices made regarding the methods described in the rest of the article.

B. Classes of methods based on input range

Reducing the range of the inputs is used in the methods described throughout the article not just to speed up their convergence, but more critically in order to ensure convergence. Methods initially designed for FP representations require the input to be in the interval [START_REF] Langhammer | Floating-Point DSP Block Architecture for FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machines[END_REF]. This is a requirement that no longer holds for this study, due to the fact that the user is allowed to specify the FxP format of the inputs.

The methods can, therefore, be devised into two separate categories. The first are methods with reduced range inputs, with an implicit assumption that the inputs are in [START_REF] Langhammer | Floating-Point DSP Block Architecture for FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machines[END_REF]. This category of methods requires a scaling of the input, in order to work with all the values of the user specified format. The second category is composed of methods that accept the full range of values of the inputs. They do not require scaling and can be used as is.

With only one exception (the case of the multipartite method), the same procedure can be used for scaling the inputs (and re-scaling the outputs at the end), for all the reduced range methods. Details of this procedure are given in what follows.

The range reduction is a classical one and consists of performing the substitution:

x = 2 s • x
where x is the input to the algorithm and x ∈ [1, 2). The value of s is chosen in such a way so that the constraints on x are satisfied. The methods presented in the next sections will take as input x .

For the post-processing phase of the range reduction, a similar substitution is applied, which results in the following set of equations:

   y = 2 -s • y y = 2 -s 2 • y y = 2 s 2 • y
where y is the output of the algorithm and y represents 1

x , 1 √ x and √

x respectively, with the semantic of x previously defined.

Figures 1 (a) and (b) present a schematic of the preprocessing and post-processing phases of the range reduction. The dotted rectangle on the top of Figure 1(a) is needed only for computing the reciprocal, and simplifies the overall architecture. The one in the bottom part of the figure computes the shift amount for the post-processing phase (it is included in this figure as it is more relevant in this context). For the postprocessing phase, the dotted rectangle represents logic that is used to align y , and adjust the shift amount accordingly, so that the output is correctly interpreted as a 2's complement number. The block labeled 'Restore Sign' is only needed for the reciprocal, and performs opposite function of the top dotted block in Figure 1(a). When needed, restoring the sign and rounding are combined into a single block.

C. Classes of methods based on approach

While the methods listed in this classification are based on the approach they take on computing the target functions, they have also been selected with the bitwidth of the input in mind. So, the possible approaches chosen are the ones which were relevant to the targeted precisions. Even so, and while not being exhaustive, the listing here covers the most used methods.

One of the main components making up contemporary FPGAs are look-up tables (LUT). This makes table-based methods an interesting and omnipresent approach. The most straightforward method is direct tabulation. It is suitable for bitwidths up to the size of the number of inputs of the LUT elements (also referred to as logic), or of the block RAM elements. The former is usually 6, while the latter is 10, or 11. The correctly rounded results can be stored directly in the table. Overflow and underflow detection can also be encoded into the stored values, without any additional resources needed.

One drawback of direct tabulation is that the resource requirement increase is exponential in the number of bits to address the table. Several methods try to overcome this shortcoming with as little extra resources as possible. These methods include the ATA method [START_REF] Wong | Fast evaluation of the elementary functions in single precision[END_REF], iATA [START_REF] Low | A memory-efficient tables-and-additions method for accurate computation of elementary functions[END_REF], the bipartite method [START_REF] Sarma | Faithful bipartite rom reciprocal tables[END_REF], SBTM [START_REF] Schulte | Symmetric bipartite tables for accurate function approximation[END_REF], STAM [START_REF] Stine | The symmetric table addition method for accurate function approximation[END_REF] and the generalized multipartite method [START_REF] Muller | A few results on table-based methods[END_REF], [START_REF] De Dinechin | Some improvements on multipartite table methods[END_REF], [START_REF] Hsiao | Table size reduction methods for faithfully rounded lookup-table-based multiplierless function evaluation[END_REF]. More details about these methods are presented in Section IV.

As the input size increases, it becomes impractical to rely solely on tabulation for implementing the target functions, even when considering the advantages brought by the multipartite methods. Another classical approach are polynomial approximation-based methods. The polynomial itself can be obtained by either using a Taylor series to approximate the function, or by using a generic polynomial approximation (like for example a minimax polynomial). Some relevant examples for the first approach would be [START_REF] Takagi | Generating a power of an operand by a table look-up and a multiplication[END_REF], [START_REF] Ito | Efficient initial approximation for multiplicative division and square root by a multiplication with operand modification[END_REF] for smaller bitwidths (discussed in Section VI-A) or [START_REF] Ercegovac | Reciprocation, square root, inverse square root, and some elementary functions using small multipliers[END_REF] for larger bitwidths (discussed in Section VII-A).

For approximating the target functions on larger bitwidths, an alternative approach to using a polynomial approximation is to use iterative methods. The most popular of these methods is the Newton-Raphson method (Section VII-B), a zero-finding method which doubles the precision of the initial guess at each iteration. To obtain an even more precise result at each iteration, there are higher order methods, with Halley's method being the next in the series after Newton-Raphson. The generalized order n version is known as the Householder method [START_REF] Borwein | Pi and the AGM: A Study in the Analytic Number Theory and Computational Complexity[END_REF]. More details on the iterative methods are shown in Section VII-C.

IV. MULTIPARTITE METHODS

The bipartite method is a well known architecture for approximating functions using only tables and no multipliers. It approximates the function using affine segments. The bipartite method can be thought of as starting with a degree 1 Taylor series in the point A:

f (x) = f (A) + f (A)(x -A)
The input interval is split into a power of 2 number of subintervals. A is selected as an easy to compute point in each interval; A can be the top k bits of x, but a better value for the approximation is obtained if the middle of each interval is selected.

The second term is indexed by the remaining bits of x that we denote by B = x -A and the derivative of the function in A, The multiplication is avoided in the bipartite method by performing a worst approximation of the derivative. Instead of storing the slope on each subinterval, one slope is stored for a larger number of subintervals:

f (A) ≈ f (C),
where C is indexed by only a few bits of A. By having this rough approximation of the slope, the entire right hand term can be tabulated and indexed by B and C. There are several optimizations regarding symmetry in this table, and the interested reader should consult [START_REF] Stine | The symmetric table addition method for accurate function approximation[END_REF], [START_REF] De Dinechin | Some improvements on multipartite table methods[END_REF].

V. APPROXIMATING ON THE FULL RANGE

In this section we will be focusing the discussions on the reciprocal function. In the bipartite architecture, the approximation quality decreases as the derivative increases. The rough approximation of the slope, especially as the input approaches zero, where the reciprocal has the asymptotic behaviour, makes the bipartite method a bad starting point for a full range architecture.

Therefore, we start-off with a degree-1 Taylor polynomial for the full-range approximation architecture. As expected, our experiments using MPFR have confirmed that the accuracy of the 1d Taylor is vastly superior to the bipartite method. Our proposed architecture will use a fine-tuned Taylor approximation for most of the input range, then use a tabulation technique for the range where the error exceeds the maximum admitted error, knowing that the total error is less than 1ulp.

It is well known that the accuracy of Taylor polynomials is significantly less than that of polynomials returned by the Remez algorithm. Open source tool Sollya [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF] provides probably the best polynomial approximations by means of the FPMinimax algorithm. We improved the accuracy of the Taylor polynomials by finding the most accurate point in a 2 dimensional space where the two dimensions were given by ulp movement of f (A) and f (A). Please note that we selected A to be the exact midpoint of each segment. The error improvement decreases the number of points which require tabulation. Since FPGA memory blocks are monolithic, they have hard thresholds. For instance lowering the number of inputs tabulated by only a few can decrease the total memory impact by half.

Depending on the input format, a large number of inputs will produce outputs which exceed the output range maximum value. For these inputs we choose to saturate the output. Additionally, a number of outputs will be in the representable range but the method error will exceed the admitted value. A tabulation method is used for these inputs. The tabulation architecture may be simple or composed of a base value and an offset (in a similar way to [START_REF] Hsiao | Table size reduction methods for faithfully rounded lookup-table-based multiplierless function evaluation[END_REF]), as depicted in Figure 2.

x 1 x Table 1 x x

Table base

Table offset Add

1 x
Fig. 2. The two architectures for the tabulation of outputs (method error exceeds the admitted value)

The resources can be lowered by sometimes tabulating the entire range: saturate range + error too large, if the total count of these inputs is less than the memory block threshold. In such case the tabulation comes for free.

The table output is equal to the output width of the architecture. If the number of table elements is large, the number of memory blocks required for implementing this table can increase significantly. Since this table stores the function value for consecutive inputs, the variation between consecutive elements in the table can be quite small. Consequently, in some cases a more efficient architecture can be obtained by sampling the input table and storing the offset values from the sample values in the offset table. Using a simple addition we can exactly reconstruct the content of the initial table, and can decrease memory impact significantly. Determining if the base+offset architecture is more efficient than simple tabulation depends on the input/output format and target FPGA memory block architecture. For one such format a designexploration stage assesses the memory impact of several decompositions: base+offset. The best architecture overall is selected by first selecting the best candidate among the base+offset architecture and the simple table-only architecture.

The output of the architecture may also underflow for a large number of inputs, also depending on the input/output format. One example is the input/output format with 16-bits of precision (unsigned) and 4 bits of fraction. For this format a separate underflow architecture is used.

The underflow architecture will use tabulation for values which are neither overflow or underflow. It has two implementation options: 1/ tabulation is used for all values which do not underflow (return 0) and 2/ tabulation is used for values which are different to 1ulp or 0. A design-space exploration phase determines the input indices for which values larger will return 1ulp (idx1ul p), and for which it would return 0 (idx0ul p). Based on these 2 values the possible architectures are explored and the one requiring the fewest memory blocks is selected. If the architectures return the same number of memory blocks, then 1/ is selected since the logic implementation is less. Figure 3 depicts the two architectures of used underflow heavy formats.

x 1 x Table Cmp. > 0 idx0 1 x x 1 x Table Cmp. > Cmp. > 0 1 idx0 idx1 1 x
Fig. 3. The two architectures of the underflow architecture

The results section presents the results of all these optimizations across an important number of input/output formats.

VI. TABULATE-AND-MULTIPLY METHOD A. First Order

One of the shortcomings of the multipartite method is the exponential growth of the memory resources required for tabulation, as the input width increases.

The method introduced in [START_REF] Takagi | Generating a power of an operand by a table look-up and a multiplication[END_REF] tries to reduce the memory requirements. It is a linear polynomial approximation of x P (for P = ±2 k , k integer), and is based on the Taylor series expansion.

The input x is split into two parts, x = x 1 + x 2 , with 1 ≤ x 1 < 2 and 0 ≤ x 2 < 2 -m . The function x P is then approximated around the point x = x 1 + 2 -m-1 , with the first two terms of the series:

x P ≈ (x 1 + 2 -m-1 ) P + (x 2 -2 -m-1 ) • P • (x 1 + 2 -m-1 ) P-1
Taking (x 1 + 2 -m-1 ) P-1 as a factor we obtain:

x P ≈ (x 1 + 2 -m-1 ) P-1 • (x 1 + 2 -m-1 + P • (x 2 -2 -m-1 )) (1)
Keeping the notations of [START_REF] Takagi | Generating a power of an operand by a table look-up and a multiplication[END_REF], we denote

C = (x 1 + 2 -m-1 ) P-1 and x = x 1 + 2 -m-1 + P • x 2 -P • 2 -m-1
The first term, C, can be obtained from a table addressed by x 1 . As for x , it can be obtained by using bit manipulations, or relatively simple operations, on x 2 and x 1 . Thus, equation ( 1) can be re-written as:

x P ≈ C • x (2)
and requires a multiplication, aside from the operations already mentioned.

In terms of the accuracy of the approximation, the error is due to the omission of the higher order terms in the Taylor series, and can be roughly approximated as:

ε method ≈ P(P -1)(x 1 + 2 -m-1 ) P-2 • 1 2 (x 2 -2 -m-1 ) 2 (3)
Which means that ε method ≈ 2 -2m-3+log |P(P-1)| . The method error ε method can be slightly improved by replacing C in eq.

(2) by C = C + P(P -1)

• x P-3 1 • 2 -2m-4
, thus eq. ( 2) becomes:

x P ≈ C • x (4)
The second term in the definition of C try to compensate for the excluded terms in the Taylor series, as shown in [START_REF] Takagi | Generating a power of an operand by a table look-up and a multiplication[END_REF].

A possible architecture for the method presented in this section is shown in Figure 4.

For the reciprocal, square root and reciprocal square root the 'x Modify' block consists mainly of inverting the bits of x 2 and concatenating them in a given order.

In terms of the implementation, the C table requires a 2 m × (w + 1)-bit table (with w the input bitwidth and m ≈ w 2 ). The multiplication is of size (w + 1) × w, and a truncated multiplier can be used so that the final result is obtained on w + 1 bits. Obtaining x requires ≈ wm LUTs, and a w-bit addition is required for the final rounding. This all makes for a quite compact architecture.

B. Second Order

The method presented in Section VI-A is based on the Taylor series, but only utilizes the first two terms of the development. As mentioned in the paragraphs discussing the method's error, the main source for the errors are the terms that are left out of the development. In what follows an extension of the method of Section VI-A is introduced, which makes use of an extra term of the Taylor series, in order to improve the accuracy of the estimation. The input is split in the same way as before, x = x 1 + x 2 , with 1 ≤ x 1 < 2 and 0 ≤ x 2 < 2 -m . The function x P is again approximated around the point x = x 1 + 2 -m-1 :

x P ≈ (x 1 + 2 -m-1 ) P + (x 2 -2 -m-1 ) • P • (x 1 + 2 -m-1 ) P-1 + 1 2 (x 2 -2 -m-1 ) 2 • P • (P -1) • (x 1 + 2 -m-1 ) P-2
Taking again (x 1 + 2 -m-1 ) P-1 as a factor we obtain:

x P ≈ (x 1 + 2 -m-1 ) P-2 • [(x 1 + 2 -m-1 ) 2 + P(x 2 -2 -m-1 )(x 1 + 2 -m-1 ) + (x 2 -2 -m-1 ) 2 P(P -2) 2 ]
For the sake of brevity and clarity, the rest of the section will deal with the case of the reciprocal (i.e. P = -1). The equations for the sqrt (P = 1 2 ) and reciprocal sqrt (P = -1 2 ) can be derived in a similar manner, and only incur different alignments for the terms, or differences in the terms stored in the tables. The previous equation thus becomes:

x -1 ≈ (x 1 + 2 -m-1 ) -3 • [(x 1 + 2 -m-1 ) 2 + (x 2 -2 -m-1 )(x 1 + 2 -m-1 ) + (x 2 -2 -m-1 ) 2 ] (5)
By regrouping the terms, eq. ( 5) can be written equivalently as:

x -1 ≈ (x 1 + 2 -m-1 ) -3 • [(x 1 + 2 -m -x 2 ) 2 + +(x 2 -2 -m-1 )(x 1 + 2 -m-1 )] (6)
Using the same style of notations as in [START_REF] Takagi | Generating a power of an operand by a table look-up and a multiplication[END_REF], we denote

D = (x 1 + 2 -m-1 ) -3 x = x 1 + 2 -m -x 2 E = x 1 + 2 -m-1 F = x 2 -2 -m-1
And eq. ( 6) becomes:

x -1 ≈ D • [(x ) 2 + E • F] (7)
The D term can be read from a table addressed by x 1 . The term E can be obtained by simply concatenating a bit of value 1 after the least significant bit (LSB) of x 1 , so it does not incur any cost in hardware. The term F can be obtained by flipping the value of the most significant bit (MSB) of x 2 , and storing the new sign of F. The term x can also be obtained by using bit manipulations, in this case it requires inverting the bits of x 2 and concatenating them to x 1 . One squaring and two additions are also required, except for the operations already mentioned. Depending on the trade-offs that the user is willing to make for the operations (or on the available hardware resources), equation ( 6) can be re-written as:

x -1 ≈ (x 1 + 2 -m-1 ) -3 • [(x 1 + 2 -m-1 ) 2 -(x 2 -2 -m-1 ) • (x 1 + 2 -m -x 2 )] (8)
Using the notations:

G = (x 1 + 2 -m-1 ) 2
Equation ( 8) becomes:

x -1 ≈ D • [G -F • x ] (9)
Compared to eq. ( 7), the method based on eq. ( 9) only requires two multiplications. On the other hand, it requires storing two coefficients in tables, as opposed to a single one. The terms D and G can be stored in a table addressed by x 1 , while F and x can be obtained by bit manipulations, as mentioned in the previous paragraphs.

The error of the method presented in this section (when computations are performed with infinite accuracy) is, as mentioned in the beginning of the section, due to the ignored terms in the Taylor series development of x P . Thus, it can be roughly approximated as:

ε method ≈ P(P -1)(P -2)(x 1 + 2 -m-1 ) -4 • 1 6 (x 2 -2 -m-1 ) 3 (10)
Which means that ε method ≈ 2 -3m-5+log |P(P-1)(P-2)/6| (with values of P of -1, 1 2 and -1 2 for the reciprocal, sqrt and reciprocal sqrt respectively). Compared to the method of Section VI-A which offers a result that is correct to approximately 2m bits, this method offers a result that is correct to approximately 3m bits. A possible architecture for implementing this method, based on equation ( 7) is shown in Figure 5. The dotted rectangle represents the part of the architecture which is different from equation [START_REF] Muller | A few results on table-based methods[END_REF]. Figure 6 presents the equivalent part for the dotted rectangle of Figure 5, in order to have a working architecture based on equation (9). Operator Modify Multiply

G Table Subtract F x / w + g / m + 1 / w + 1 / 2m + 2
Fig. 6. Alternate Second Order Tabulate-and-Multiply Method Architecture

In terms of the implementation, the squaring of x and the final multiplication need not be performed on the whole precision of their inputs. Truncating these operations, on the other hand, entails rounding errors, which may prevent from reaching the target output precision (faithfully rounded result). So, a larger intermediary precision is used, adding a number of g extra guard bits to the datapaths.

The total error budget (ε total ) is divided as follows: ε total = ε f inal round + ε method + ε round , with a final goal of ε total < 2 -w . Half of an ulp (unit in the las place) is allocated to the final rounding (ε f inal round ). The rest is divided between the method error (ε method ) and the rounding errors (ε round ), which are due to the approximations performed during the computations. Determining ε method will give the value that needs to be used for m, while determining ε round will determine the value of g and the widths of the datapaths.

In order to determine ε round , eq. ( 7) is re-written to reflect the approximations:

x -1 ≈ D • [ (x ) 2 + E • F]
and can be further re-written as:

x -1 ≈ (D + ε table ) • [(x ) 2 + ε sqr + E • F]
where the tilde terms are the approximations that are obtained from the table and the squarer, respectively, and the ε are the corresponding errors. Developing the equation results in:

x -1 ≈ D • [(x ) 2 + E • F] + ε table [(x ) 2 + F • E] + ε sqr (D + ε table ) + ε mult
from which ε round can be identified as:

ε round = ε table [(x ) 2 + F • E] + ε sqr (D + ε table ) + ε mult
Using the upper bounds for the expressions (x ) + E • F and D results in:

ε round ≤ 4ε table + ε sqr + ε sqr ε table + ε mult (11)
If the same internal precision is used for all the operations, meaning that ε table = ε mult = ε sqr = 2 -w-g , and keeping in mind that ε round < 2 -w-2 , the following expression is found for g:

g > 2 + log(6 + 2 -w-g )
This is, however, a pessimistic value, and testing has shown that lower values of g might also suffice. It is also confirmed by considering individual precisions for the datapaths (given by the ε terms of eq. ( 11)). If eq. ( 9) is used for implementing the method, a similar approach can be used in order to determine ε . As previously, eq. ( 9) is then re-written as:

x -1 ≈ D • [G -F • x ]
and then as:

x -1 ≈ (D + ε table ) • [G -F • x + ε mult1 ] + ε mult2
and finally as:

x -1 ≈ D • [G -F • x ] + G • ε table -F • x • ε table + D • ε mult1 + ε mult1 • ε table + ε mult2
from which ε round can be identified as:

ε round = G • ε table -F • x • ε table + D • ε mult1 + ε mult1 • ε table + ε mult2
Using the same reasoning and the same constraints as before, a relation is found for g, as well:

g = 2 + log(4 -2-m -1 + 2 -w-g )
The same remarks can be made for the slightly pessimistic relations found for g here, as the ones made for the previous architecture, in the previous paragraphs. The methods presented in Section VI can tackle, from a practical point of view, precisions of up to about 24 bits.

VII. LARGE PRECISION METHODS

A. Taylor-based Method

The methods introduced in Section VI rely on first and second order approximations using the Taylor series (as a side note, the multipartite methods could be seen as well as first order polynomial approximations). As the input width increases, the methods require either an increasing number of terms, as well as a rapidly increasing amount of storage for the coefficients.

The method presented in [START_REF] Ercegovac | Reciprocation, square root, inverse square root, and some elementary functions using small multipliers[END_REF] tries to overcome these obstacles with an approach consisting of three steps. The first step consists of a range reduction of the input. The next step is the evaluation of the function on the reduced argument using a series expansion of the target function. Last is a postprocessing phase, which is required due to the first step.

The main premise of the method is a classical one: the evaluation of the function on the reduced argument is less expensive than on the original input. This is due to the higher order terms in the evaluation of the series being shifted out to lower weights than the target precision, thus they do not need to be computed.

The range reduction transforms the input 1

≤ x < 2 into -2 -k < x < 2 -k . x is obtained as x = x • x -1, where x is a k + 1 bit approximation of 1
x . x can be obtained from a table addressed by the k MSBs of x.

The evaluation step consists of approximating the target function ( f ) at x . f can be expressed as:

f (x ) = C 0 +C 1 • x +C 2 • (x ) 2 +C 3 • (x ) 3 + . . .
where the coefficients C i are given by the Taylor series development. The input x is split into 4 chunks of size k, x 1 to x 4 ; as |x | < 2 -k , x 1 can be left out, as it only contains sign bits. Replacing x in the previous equation and eliminating terms of weight of less than 2 -4k results in:

f (x ) ≈ C 0 +C 1 • x +C 2 • x 1 • 2 -4k + 2C 2 • x 2 • x 3 • 2 -5k +C 3 • (x 2 ) 3 2 -6k
The post processing consists of multiplying the result obtained at the previous step with a term that can be obtained similarly as x. This term is x, 1 √ x and

√

x for the reciprocal, sqrt and reciprocal sqrt.

One of the advantages of this method would be the small multipliers and tables required for implementation, which seem to map quite well to the hardware resources available on the current generation FPGAs. Another advantage would be the simple range reduction and post-processing phase.

It would be interesting to see, however, to what degree the method could be improved by replacing the approximation of f (x ) by the Taylor series development with a more accurate polynomial approximation.

B. Newton-Raphson Method

Probably the best known iterative method is the Newton-Raphson scheme. This is a root-finding scheme and has quadratic convergence. In order to deduce the method, the starting point is again the Taylor series, around the point x n , limited to the first two terms:

f (x) = f (x n ) + f (x n ) • (x -x n ) A root of f (x) satisfies f (x) = 0, so: f (x n ) + f (x n ) • (x -x n ) ≈ 0 which gives: x ≈ x n - f (x n ) f (x n )
Thus we can obtain the recurrence relation:

x n+1 = x n -f (x n ) f (x n ) .
The error entailed by this iteration can be approximated by a recurrence relation as well:

ε n+1 ≈ ε 2 n • f (x n ) 2 • f (x n )
which shows that the method converges quadratically.

Again, for the sake of brevity, the following discussions will be limited to the case of the reciprocal. In order to compute the reciprocal of a number a, a function f is needed so that f

( 1 a ) = 0. A convenient choice for the function f is f (x) = 1
xa. Replacing f into the recurrence relation results in the following scheme:

x n+1 = x n • (2 -a • x n ) (12)
In terms of the implementation, the iteration given by eq. ( 12) only requires additions and multiplications, once an initial approximation has been obtained.

We perform an error analysis similar to those of Section VI-B. The total error is again divided as ε total = ε f inal round + ε round + ε method , and ε total < 2 -w . Eq. ( 12) thus becomes:

x n+1 = x n • (2 -a • x n ) (13)
This assumes that we have an initial approximation correctly rounded to m = w 2 bits (which implies an error less than 2 -m-1 ). Writing explicitly the errors in eq. ( 13), it becomes:

x n+1 = (x n + 2 -m-1 ) • (2 -a • (x n + 2 -m-1 ) + ε mult1 ) + ε mult2
which can be written as:

x n+1 = x n • (2 -a • x n ) + ε round
where ε round represents the rounding errors

ε round = 2 -m-1 • (2 -a • x n ) -2 -m-1 • a • (x n + 2 -m-1 ) + ε mult1 • (x n + 2 -m-1 ) + ε mult2 (14)
The sum of the first two terms of eq. ( 14) is of the order of magnitude of 2 -2m-2 . If ε mult1 and ε mult2 are of the form 2 -w-g (where g gives the number of extra guard bits), eq. ( 14) shows that at least 3 guard bits are needed. A good option for the initial approximation is the use of a bipartite/multipartite table. This is also the approach used in this article. For target precisions of up to approx. 32 bits w 2 falls in the range of values that are suitable for use with the bipartite method. Increasing the number of iterations is also possible, thus increasing the number of multiplications, but decreasing of resources needed to obtain the initial approximation.

C. Generalized Iterative Methods

The Newton-Raphson method discussed in Section VII-B has quadratic convergence. However, there are methods with faster convergence. The second order Newton-Raphson method is also known as the Halley method. In order to deduce the method, we will have the same approach as in the previous section, starting with the Taylor series around the point x n , but using the first three terms:

f (x) = f (x n ) + f (x n ) • (x -x n ) + 1 2 f (x n ) • (x -x n ) 2
A root of f (x) satisfies f (x) = 0, so:

f (x n ) + f (x n ) • (x -x n ) + 1 2 f (x n ) • (x -x n ) 2 ≈ 0
which can be re-grouped as:

x n+1 = x n - f (x n ) f (x n ) + f (x n )•(x n+1 -x n ) 2 Knowing that x n+1 -x n = -f (x n ) f (x n )
results in the iteration formula:

x n+1 = x n - 2 • f (x n ) • f (x n ) 2 • ( f (x n )) 2 -f (x n ) • f (x n ) (15)
Equation ( 15) is commonly known as the Halley iteration/method. The error entailed by the method can be expressed a recurrence as well, as in the case of the Newton-Raphson iteration:

ε n+1 ≈ ε 3 n • f (x n ) 6 • f (x n )
which shows the cubic convergence of the method. The Newton-Raphson and the Halley method are the first and the second in a class of iterative methods known as the Householder methods [START_REF] Borwein | Pi and the AGM: A Study in the Analytic Number Theory and Computational Complexity[END_REF]. The n + 1-st term in the class has the form: [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF] where the (n) superscript represents the n-th order derivative. It can be easily verified that choosing n = 0 results in the Newton-Raphson iteration, while choosing n = 1 results in the Halley iteration. In terms of the error, the Householder method progresses at a pace of order n+1 towards the correct solution.

x n+1 = x n + (n + 1) • 1 f (x n ) (n) 1 f (x n ) (n+1)
The Halley method (or the Householder method, more generally) replaces the tangent to the function plot by a curve (or a higher order curve, in the case of Householder). This curve has a higher number of derivatives in common with the function plot, at the point of the approximation. This should, in principle, better fit the plot, giving, thus, a better approximation. However, as remarked in [START_REF] Flynn | On division by functional iteration[END_REF], in the case of the reciprocal this expansion is not particularly useful. Plugging-in the same function as in the case of the Newton-Raphson iteration in the Halley iteration requires the computation of the inverse that we are trying to approximate in the first place.

On the other hand, as remarked in [START_REF] Rabinowitz | Multiple-precision division[END_REF] and even further back presented in [START_REF] Willers | Practical analysis: graphical and numerical methods[END_REF] taking a different approach to obtaining the iteration in the first place seems to be more effective.

The starting point for the alternative approach is again the Taylor series around the point x n :

f (x) = f (x n ) + f (x n ) 1! • (x -x n ) + f (x n ) 2! • (x -x n ) 2 + f (x n ) 3! • (x -x n ) 3 + . . . (17)
Again, a root of f (x) satisfies f (x) = 0. The value of xx n is expressed as a power series of f (x n ):

(x -x n ) = a • f (x n ) + b • ( f (x n )) 2 + c • ( f (x n )) 3 + . . . (18)
By replacing eq. ( 18) in ( 17) and using the fact that f (x) = 0 results in:

0 = f (x n ) + f (x n ) 1! • (a • f (x n ) + b • ( f (x n )) 2 + c • ( f (x n )) 3 + . . .) + f (x n ) 2! • (a • f (x n ) + b • ( f (x n )) 2 + c • ( f (x n )) 3 + . . .) 2 + . . . (19)
If eq. ( 19) is seen as an equation with f (x n ) as a variable, the coefficients of the same powers of f (x n ) on the two sides of the equation can be identified. Thus, the a, b, c, . . . coefficients can be found:

a = - 1 f (x n ) b = - f (x n ) 2 • ( f (x n )) 3 c = - ( f (x n )) 2 2 • ( f (x n )) 5 . . . (20)
Replacing the values of eq. ( 20) in eq. ( 18) results in an iteration of the form (showing only the first three terms):

x n+1 = x n - f (x n ) f (x n ) - ( f (x n )) 2 • f (x n ) 2( f (x n )) 3 -. . . (21)
Eq. ( 21) can be used as an alternative to Halley's method.

In order to obtain the cubic iteration for the reciprocal function, f (x) in eq. ( 21) is chosen as f (x) = 1

xa. This results in:

x n+1 = x n • (1 + h n (1 + h n )) h n = 1 -a • x n (22)
In order to obtain the iteration for the reciprocal sqrt, f (x) in eq. ( 21) is chosen as f (x) = 1

x 2a. This results in:

x n+1 = 1 8 • x n • (8 + h n (4 + 3 • h n )) h n = 1 -a • x 2 n
Both the iteration for the reciprocal and the one for reciprocal sqrt have cubic convergence.

The rest of this section is dedicated to the error analysis for the datapath implementing eq. ( 22). The evaluation, on the other hand, is done as (where h n keeps the meaning of eq. ( 22)):

x n+1 = x n + x n • (h n + h 2 n ) (23)
This form of the iteration takes advantage of the fact that h n < 2 -m (the proof and reasoning are similar to the range reduction of [START_REF] Ercegovac | Reciprocation, square root, inverse square root, and some elementary functions using small multipliers[END_REF]). This means that less bits are needed for the squaring of h n and for the multiplication x n • (h n + h 2 n ), due to the terms being shifted to the right by 2m and m bits respectively.

As for the case of the previous error analyses, the total error is divided in three parts, ε total = ε f inal round + ε round + ε method < 2 -w . The final result is rounded, so ε f inal round < 2 -w-1 . The method error, ε method , was previously discussed. What remains to be determined are the rounding errors, ε round . These are due to the initial approximation and to truncating the multiplications in eq. ( 23). Making these errors explicit, eq. ( 23) becomes (where m is the precision of the x n approximation):

x n+1 = (x n + 2 -m-1 ) + (x n + 2 -m-1 ) • (h n -a • 2 -m-1 -ε mult + (h n -a • 2 -m-1 -ε mult ) 2 ) or x n+1 = x n + x n • (h n + h 2 n ) + ε round where ε round = 2 -m-1 • (1 + h n + h 2 n )+ + (x n + 2 -m-1 ) • (a • 2 -m-1 + ε mult )• • (a • 2 -m-1 + ε mult -1 -2 • h n )
which can be re-written as:

ε round = 2 -m-1 • (1 + h n + h 2 n )+ + (x n + 2 -m-1 ) • [((-a)2 -m-1 + a 2 2 -2m-2 -a2 -m h n )+ + ε mult (1 -ε mult -2h n + a2 -m )] (24)
where ε mult denotes the rounding error due to a multiplication or of a squaring. The terms of equation (24) that do not contain ε mult will result in a term that is of the order of 2 -3m-3 . Thus, if ε mult is of the form 2 -w-g , where g is the number of guard bits, we should ensure that g > 3 for a faithfully rounded result. This has also been confirmed through testing.

VIII. RESULTS

In this section we present the resource utilization and performance of cores generated using the previously described methods.

In Table II we present the resources and performance of the extended low order Taylor method that we introduced in Section V for varying input formats. The table is split into two sides: the left side presents the performance results for the default Taylor series coefficients whereas the right side presents the results for enhanced coefficients, found starting from the previous coefficients. The search methods applied allows exploring the design space and returning the architecture which approaches the functions best in less than 1 second, on recent compute architecture (Core i7).

There are several architectural variations of this solution:

• specialized underflow architecture; this architecture will generally trigger when the formats are of such sort that a significant number of inputs will underflow. There are 2 variations of this architecture: as depicted in Figure 3. An example which triggers the left architecture is the [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF][START_REF] Wong | Fast evaluation of the elementary functions in single precision[END_REF] format, and which triggers the right is the (16,5) format.

• architecture for which the approximation is sufficiently accurate on the compute range (except overflow). An example of this architecture is given by the [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF][START_REF] Borwein | Pi and the AGM: A Study in the Analytic Number Theory and Computational Complexity[END_REF] format.

• architecture for which an extra table is used for inputs for which the 1st degree approximation is not sufficiently accurate. The architecture of the extra table is selected in order to minimize memory block utilization. This table may also capture the inputs which cause the output to overflow, if doing so does not increase the memory block requirements. For the table itself there are two subarchitectures possible, depicted in Figure 2. For instance, format [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF][START_REF] Takagi | Generating a power of an operand by a table look-up and a multiplication[END_REF] will save 1 memory block by using the base+offset architecture, and will have the table indexed by x. Format [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF][START_REF] Ercegovac | Reciprocation, square root, inverse square root, and some elementary functions using small multipliers[END_REF] will also save 1 memory block but the table will be indexed by an offset version of x; the range close to 0 will be handled separately as handling it in the table would increase the memory requirements.

Among the presented formats, some show better savings than others in therms of the optimizations effect. For the ones which show little-to-no savings, we observed this was due to the size of memories; we were unable to reduce the number of data elements sufficiently (usually below a power of 2) to trigger a smaller memory block utilization.

The final row in Table II shows the resource requirements and performance for a bipartite implementation which uses an input shifter and leading-zero-counter to normalize the input, and a right shifter to denormalize the result. This type of implementation is regarded as classical when dealing with fixed-point implementations on a generic range. This implementation is generally agnostic to the input/output format, so we only give results for one such format. When comparing this method with the previously optimized method for the whole range we observe that the latency of this implementation is significantly larger. Next, since based on the bipartite implementation, the architecture does not consume any DSP blocks but it consumes more logic (ALMs) than all the presented architectures. In terms of memory blocks, the architecture performs significantly better for inputs ranging from [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF][START_REF] De Dinechin | Some improvements on multipartite table methods[END_REF] to [START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF][START_REF] Ercegovac | Reciprocation, square root, inverse square root, and some elementary functions using small multipliers[END_REF] but is on par, or consumes more memories for the rest of the formats.

These results show that in the context of a core generator, the user requirements will determine which architecture will be returned. If the objective is latency or ALMs alone then the method for the entire range will be better; if targeting memory blocks, then the tie-break will be played on formats.

In Table III we present the resource requirements together with the performance for cores having inputs in the interval [START_REF] Langhammer | Floating-Point DSP Block Architecture for FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machines[END_REF], and precision varying from 10 bits up to 32 bits. The varying precisions will trigger the different architectures: table-based, bipartite, Taylor-based, and high-order based on Newton-Raphson, Halley and Taylor. As expected, bipartite methods perform well up to 16-bit precisions. The method Tabulate and Multiply presented in Section VI-A performs well up to 24 bits. With higher precision one can observe the increase in memory blocks for this architecture. Finally, for 32-bits we compare the requirements and performance of a Newton-Raphson-based implementation bootstrapped by a bipartite approximation, a higher-order Taylor based implementation, and two variations of Halley's method: first bootstrapped by a table, and the second bootstrapped by a bipartite approximation. These implementations show interesting tradeoffs: Newton-Raphson+ bipartite shows the highest number of memory blocks (5) but the lowest number of DSPs (2), an average-low number of ALMs, and a short latency; Taylor uses the most ALMs, but fewest memory blocks, highest number of DSPs and highest latency; Halley bootstrapped by a table has the lowest ALM and latency but highest DSP and memory; Halley bootstrapped by bipartite reduces the number of memory blocks to minimum, but slightly increases the number of ALMs.

Table III misses the implementation of the second-order Tabulate and Multiply method presented in section VI-B which is expected to lower memory requirements around the 24-bit precisions where the current first order implementation uses 5 memory blocks. The results for this method will be present in the final version of the paper.

Table IV presents the resource utilization and performance of small precision floating-point cores. The utilization numbers are presented in order to contrast these implementations against the fixed-point counterparts. The floating-point implementations are based on a piecewise-polynomial approximation with half-precision using a degree 1 polynomial, whereas remaining precisions presented in the table using a degree 2 polynomial. Implementations use the Horner scheme and the coefficient sizes follow closely [START_REF] De Dinechin | Automatic generation of polynomial-based hardware architectures for function evaluation[END_REF]. Horner scheme also uses truncated multipliers, however the low degree polynomials and larger monolithic multipliers in the target FPGA device diminish the benefits. Please note that the kernels used in these implementations directly use normalized inputs, so are comparable to Table III.

IX. CONCLUSION

In this article we have revisited some well known methods for implementation of fixed-point kernels, and depicted their performance on a contemporary FPGA device: an Altera

Fig. 1 .

 1 Fig. 1. Scaling range reduction

TABLE I SUBSET

 I OF DSP-BLOCK FEATURES FOR CYCLONEV DEVICES

	Blocks Mode	CycloneV
		9 × 9	3
		16 × 16	2
	1	18 × 18 → 32 18 × 18	2 2
		18 × 19	2
		27 × 27	1
		2 18 × 18 MADD	1
		2 18 × 19 MADD	1
		18 × 18 + 36-bit	1
	2	2 27 × 27 MADD	1

  1.xxx . . . . . . x xxx . . . . . . x

	x 1	x 2
	/ m	x
	C Table	x Modify
		/ w + 1
		/ w + g
		Multiply
		/ w + g
		Round
		x P
	Fig. 4. Tabulate-and-Multiply Method Architecture

  1.xxx . . . . . . x xxx . . . . . . x

	x 1		x 2
	/ m	x
	D Table		x Modify
		E	F	x
	/ w + g	Multiply	SQR
				/ w + g
			Add
			/ w + g
		Multiply
		/ w + g
		Round
		x P

/ m + 1/ m + 1 / w + 1 / 2m + 2

Fig. 5. Second Order Tabulate-and-Multiply Method Architecture

ACKNOWLEDGEMENTS

This work was supported by the Altera European Technology Centre (that Matei Is ¸toan was part of at the moment when the works described in this paper started) and the French National Research Agency through the INS program MetaLibm.

CycloneV. Unlike previous works which generally focus implementations of normalized inputs ([1,2) or [0.5,1)), we have focused in this work on the the full implementation range.

Where possible we have modified the classical architecture to directly handle the full input range; the results of this work were competitive architectures which trade off logic resources for memory blocks. Our results also show that the behaviour and performance of the modified full-range architectures are highly dependent on the input and output format. Therefore, our architectures are more likely to be regarded as implementation tradeoffs. We have shown our results on the reciprocal function, but we have used the same templates for the reciprocal square root and square root functions. We have also revisited higher order methods based on the Newton-Raphson algorithm, Taylor series, and Halley's method, and proposed combined architectures; there, Newton-Raphson and Halley are bootstrapped using a bipartite-based implementation. The results shown in this paper confirm that fixed-point arithmetic is very relevant for contemporary FPGAs containing up to thousands of embedded memory blocks and multipliers, and the low resources of these cores in combination with the low resource of primitive components will keep fixed-point datapaths relevant in resource critical scenarios, and when expertise exists to bound variable ranges.