Square root in fixed point VHDL

Square root in fixed point VHDL

The most common math algorithms that are needed when we deal with numbers in embedded systems are sum/subtraction, multiplication, division. When we move from scalar to vector quantities and complex numbers, we invariably run into square root when we operate with lengths of vectors. Square root also comes up as the geometric mean of numbers as well as different energy and rms calculations of measurements. Hence it is next most important as far as elementary math operations are concerned

In this post we will develop the basic algorithm for parametrizable word length and accuracy of a fixed point square root using the newton Raphson inverse square root algorithm. This algorithm will then be tested against the real_math vhdl package implementation, simulated in ghdl and finally synthesized to Efinix Titanium FPGA. The efinix titanium sources are available in the example projects github page and the multiplier and square root fixed point algorithms shown here are available in the hVHDL fixed point library . All of the testbenches in the project are handled by VUnit and the simulations can be launched with

python vuint_run.py -p 32 --gtkwave-fmt ghw The waveforms can then be found in the vunit_out/test_output/test_name /ghdl/wave.ghw

High level design of vhdl modules

Although VHDL is used to describe functionality of digital circuits it has several high level language features like custom and composite types and overloading of subroutines. With these features we can create interfaces to our code which allow us to make our modules easier to use. Since digital circuits are massively paralleled and since we are implementing fixed point algorithm the important points of abstraction are timing and radix handling. Making subroutines that handle these allows us to concentrate on using the algorithms instead of worrying about the low level details.

The following snippet shows how we use recods and subroutines to create interfaces for VHDL modules. The create_multiplier and create_square_root procedures set up logic around the registers that are defined in the records and subsequent request_square_root procedure then modifies the state to trigger the calculation. Notice that there is no counters outside the square root, instead the record-procedure pair that contains all of the state machines and timing inside the square root itself and we read the instant when calculation is ready by polling square_root_is_ready function.

-- architecture of a module signal multiplier : multiplier_record := init_multiplier signal square_root : square_root_record := init_square_root; begin process_with_fixed_point_square_root : process(clock) begin if rising_edge(clock) then create_multiplier(multiplier); create_square_root(square_root, multiplier); if square_root_calculation_is_needed then request_square_root(square_root, input_number); end if; if square_root_is_ready(square_root) then result <= get_square_root_result(square_root, radix); end if; end process; The interface of the sqrt module has and interface with handshakes.The create_procedure creates the module or initializes it, then the calculation is launched by request_square_root procedure and lastly we wait for the module to be ready and then get_square_root_result and with some desired radix.

Using interfaces to modules from which the result can be fetched with desired radix, we substantially simplify the fixed point calculations.With the interface we don't track the radix outside of the code, but instead we choose what radix we use and the module returns its output with the desired radix.With the square_root_is_ready function, we hide the timing information from the application code as the module it self keeps track of it own timing.This helps to contain changes inside the interface as if we use this interface and change the latency for example to add registers to meet timing, the application code will adapt automatically.--> <!--

We create an interface using records and subroutines, and the square root functionality is implemented without the need for a separate entity; instead, we use records and subroutines directly. If desired, we could encapsulate it within an entity, allowing us to manage input and output data through separate record structures.he code snippet below shows how we can use this high level abstraction in VHDL modules. We write an interface with records and subroutines and the square root here is designed without a entity directly with records and subroutines. If we put it inside an entity, then we would have records for the input and output data separately. --> <!--

As seen in the snippet square root functionality will be implemented with records and subroutines directly without a bounding entity. This allows the square root calculation to be added into any process where it is needed with only a few lines of code and it can be tested and designed with minimal dependencies. If desired, we could encapsulate the square root module within an entity, in which case we would modify the interface to use input and output data through separate records.

Algorithms for square root

For a fast and accurate implementation of a square root there are two similar algorithms, the Goldschmidt algorithm and inverse root finding algorithm based on Newton Raphson iteration. An article on Goldschmidt Algorithm , by Michael Morris has been posted previously on embeddedrelated and the fast calculation of inverse square root with NR is legendary enough to have its own wikipedia article . Both of these algorithms are iterative which means that we can just run run them in a loop by using the result of previous run to get more accurate result with the next one. We will quickly take a look at the similarities of these two algorithms.

Goldschmidth square root

Goldschmidt estimates both the $\sqrt{x}$ and $\frac{1}{\sqrt{x}}$ and the algorithm allows parallel or pipelined calculation of two of the multiplications, and it requires three initial values to get started.

$$x_{0} = S\left(\frac{1}{\sqrt{S}}\right)_{est}$$

$$h_{0} = \frac{1}{2}\left(\frac{1}{\sqrt{S}}\right)_{est}$$

$$r_{0} = \frac{1-2x_0 h_0}{2}$$

The actual algorithm which converges to both $\sqrt{S} = x_i$ and $\frac{1}{\sqrt{S}} = 2h_i$ is

$$x_{i} = x_{i-1}\cdot (1+r_{i-1})$$

$$h_{i} = h_{i-1}\cdot (1+r_{i-1})$$

$$r_{i} = \frac{1-2x_i \cdot h_i}{2}$$

Calculation of $x_i$ and $h_i$ can be done in parallel and their result is needed to calculate $r_i$. The algorithms have the mandary multiplications written out to show where we must use an actual multiplier.

Newton-Raphson inverse square root

Using the NR root finding algorithm $x_{n+1} = x_n - \frac{f(x_n)}{f'(x)}$  to $\frac{1}{x^2}-S = 0$ which has its root at $\frac{1}{\sqrt{x}}$  the gives the classic recursive NR square root algorithm

$$x_{n+1} = 1.5x_n - 0.5Sx_n^3$$

With a bit of rewriting we can see the similarities between the Newton-Raphson and Goldschmidt algorithms

$$S_{xi} = S\cdot x_{i-1}$$

$$x_{xi} = x_{i-1} \cdot x_{i-1}$$

$$x_{i} = 1.5x_{i-1} - 0.5x_{xi}\cdot S_{xi}$$

To get the NR iteration started we only require an estimate of the inverse square root of $S$

$$x_{0} = \left(\frac{1}{\sqrt{S}}\right)_{est}$$

The NR algorithm gives the inverse square root from which we can easily get $\sqrt{x}$ as $x \cdot \frac{1}{\sqrt{x}} = \sqrt{x}$ . We can see from the rewritten equations the two algorithms have a very similar structure. Both have two pipelineable multiplications and require a third which multiplies the result of the previous two. Due to the easier initialization of the algorithm we will use the NR for our square root. The hVHDL library also has a division algorithm that uses practically identical implementation.

Fixed point implementation

The algorithm presented here works equally for fixed or floating point implementation. The difference between fixed and floating point is that in fixed point we code the exponent, or equivalently the number of integer bits, into our implementation since it is not dynamically carried with the number representation as is the case with floats.

Choosing radix of NR algorithm

Choosing the number of integer bits for the implementation of the inverse square root algorithm is quite straightforward. Since the convergence is fast as long as the initial guess is close enough, we will restrict the input range to be within $[1.0, 2.0)$ and the inverse square root of this range maps to $[0.707, 1.0]$. We will use 2 integer bits for the algorithm as this is the minimal amount integer bits that cover both ranges therefore our radix will be the determined by word length - 2. We will use this radix throughout the calculation and thus we don't need to adjust any of the radixes or wordlengths during the implementation of the iteration.

With this choice, the bit vectors that represent the minimum and maximum values are $1.00000...$ and $1.111111...$. Hence we can map any bit vector that comes in by left shifting it until it has no leading zeros and post scale by left shifting again with half of the number of leading zeros in the input number.

NR iteration module

To start the implementation we can break down the NR iteration to single operations of additions and multiplications to get

step 0 : x <= initial_value; Loop: step 1 : ax <= a * x; three_halfs <= x + x/2; step 2 : x2 <= x * x; step 3 : ax3 <= ax * x2; step 4 : x <= three_halfs - ax3; goto Loop step 5 : result <= x; Steps 1 through 4 can be run in a loop for greater accuracy. We will use a fixed point multiplier module to do the multiplications. This module is similarly written in a record and a procedure and it has parametrizable input and output pipeline depths. The module is parametrized with the use of configuration packages for the bit widths and pipeline depths. We could do the parametrization with package generics however some of the tools that we still use with fpgas, mainly intel quartus and older versions of vivado, do not allow the use of VHDL2008 features.

The registers that we need for building the inverse square root are defined in a isqrt_record. The two state counter registers state_counter and state_counter2 are used for sequencing the calculations, registers for the loop_value allows us to keep track of how many times we want to run the NR iteration loop. The full inverse square root record is as follows

type isqrt_record is record x_squared : signed(int_word_length-1 downto 0); threehalfs : signed(int_word_length-1 downto 0); x : signed(int_word_length-1 downto 0); result : signed(int_word_length-1 downto 0); sign_input_value : signed(int_word_length-1 downto 0); state_counter : natural range 0 to 7; state_counter2 : natural range 0 to 7; isqrt_is_ready : boolean; loop_value : natural range 0 to 7; end record; We build the actual logic that creates the inverse square root algorithm with a create_isqrt procedure. This procedure gets the multiplier through the argument list as this allows reusing the same hw multiplier outside of the create_isqrt procedure for the rest of the square root algorithm.

procedure create_isqrt ( signal self : inout isqrt_record; signal multiplier : inout multiplier_record ) is variable mult_result : signed(int_word_length-1 downto 0); variable inverse_square_root_result : signed(int_word_length-1 downto 0); begin CASE self.state_counter is WHEN 0 => multiply_and_increment_counter(multiplier,self.state_counter, self.x, self.x); self.threehalfs <= self.x + self.x/2; WHEN 1 => multiply_and_increment_counter(multiplier,self.state_counter, self.x, self.sign_input_value); WHEN others => --do nothign end CASE; self.isqrt_is_ready <= false; CASE self.state_counter2 is WHEN 0 => if multiplier_is_ready(multiplier) then self.x_squared <= get_multiplier_result(multiplier, isqrt_radix); self.state_counter2 <= self.state_counter2 + 1; end if; WHEN 1 => if multiplier_is_ready(multiplier) then multiply(multiplier, self.x_squared, get_multiplier_result(multiplier,isqrt_radix)); self.state_counter2 <= self.state_counter2 + 1; end if; WHEN 2 => if multiplier_is_ready(multiplier) then mult_result := get_multiplier_result(multiplier,isqrt_radix+1); inverse_square_root_result := self.threehalfs - mult_result; self.x <= inverse_square_root_result; self.state_counter2 <= self.state_counter2 + 1; if self.loop_value > 0 then request_isqrt(self , self.sign_input_value , inverse_square_root_result , self.loop_value); else self.isqrt_is_ready <= true; self.result <= inverse_square_root_result; end if; end if; WHEN others => --do nothign end CASE; end create_isqrt; If our multiplier has only one pipeline cycle, then the first multiplication is ready at the same time as the second one is pipelined into the multiplier.

By having two state counters we can handle this situation.We also schedule the calculation of $1.5\cdot x$ in the first state_counter.

The isqrt calculation is launched by setting the input value and by resetting the state counters with request_isqrt procedure. Both of the state machines run through their steps in a linear sequence and halt when they run to completion. The request sets both state machines to zero state which also resets the state to known and defined position.

procedure request_isqrt ( signal self : inout isqrt_record; input_number : signed; guess : signed; number_of_loops : natural range 1 to 7 ) is begin self.sign_input_value <= input_number; self.state_counter <= 0; self.state_counter2 <= 0; self.x <= guess; self.loop_value <= number_of_loops-1; end request_isqrt; The full implementation of the inverse square root with its interface functions can be found in fixed_isqrt_pkg.vhd

The algorithm is developed in a testbench that can be found here . The algorithm converges if the initial value is close to the real root. As noted before the algorithm converges quadratically, which approximately means that we double the number of correct bits every iteration. This convergence is so fast that for input number between 1.0 and 2.0 a constant initial value is good enough to give roughly 32 bits of accuracy in just 4 iterations. The error behavior with 4 iterations is shown in Fig 1. The error is calculated between the fixed point square root and a reference value that is obtained by counting square root with the real_math library. The figure shows maximum error which with 4 iterations is roughly $1.8\cdot 10^{-10}$.

multiple initial values

Since the convergence is very fast and the algorithm only requires 3 multiplications of which two can be pipelined, we have an interesting situation with the initial value. There is no benefit obtained from getting a better initial guess unless it is faster to obtain than what it takes to run another iteration. Even a really crude guess gets us fast convergence and the higher accuracy we desire, the higher accuracy initial guess we will need to avoid an iteration through the algorithm. Fig 2. shows the error behavior with the initial value being one of 16 possibilities. This allows us to only use 2 iterations for error less than $2.0\cdot 10^{-7}$ which is close to accuracy of a 32 bit float.