How to perform small value approximation for sqrt(x) on FPGA

Question

I am trying to implement a fixed point routine which involves computing the value of \$ \sqrt{x} \$ for small \$x\$ that approaches \$0\$. The target architecture is an FPGA. One problem is that this function does not lend itself easily towards the use of Taylor's expansion. One can see that for small values of x, the slope of \$\sqrt{x}\$ goes to infinity when \$x\$ approaches \$0\$, therefore evaluating the function using a power series involves multiplying huge coefficients with a small \$x\$. This method is therefore numerically unstable.

Using an iterative approach, the Newton-Raphson yields the following iterative equation: \$x_{n+1} = \frac {x_{n}}{2}- \frac{\alpha} {2x_{n}}\$, where we are trying to approximate \$\sqrt {\alpha}\$. But once again, since \$\alpha\$ is small, \$x_{n}\$ would likewise have to be small for the solution to converge. Since the equation involves dividing a small number by another small number, chances are that fixed point arithmetic would fail.

With that, I would like to know how to implement small value approximation for \$\sqrt{x}\$ using fixed point arithmetic, either using precomputated coefficients or iterative methods.

If you are targeting an FPGA, the first and most important question is which precision you want. You say you want to use fix point: which precision for the input, which precision for the result? In fixed point (as in integers) there is no "approaching zero". There is just a smallest number you are interested in. — Philippe, Mar 14 '12 at 08:26

score 5 · Accepted Answer · answered Nov 02 '11 at 16:15

A routine which I have used before (I don't know if it's a "proper" one or not) is a divide-and-conquer approach.

You start off with an arbitrary upper and lower value (say 5 and 0 respectively - the highest and lowest square roots you want to find) and find the mid-point between them. Square that value.

If the squared value is greater than your target, set the upper value to be your squared value. If it's lower, set the lower value.

Repeat until either the square value matches your lookup value, or you have executed enough iterations to be as accurate as you like.

Here's a little version I have knocked together in perl:

#!/usr/bin/perl

my $val = shift;

my $max = 5;
my $min = 0;

my $iterations = 0;
my $maxiter = 40;

while(($max > $min) and ($iterations<$maxiter))
{
    $iterations++;
    my $diff = $min + ($max - $min) / 2;
    my $square = $diff * $diff;

    if($square == $val)
    {

        print "Square root found at $diff\n";
        print "$iterations iterations\n";
        exit(0);
    } else {
        if($square > $val)
        {
            $max = $diff;
        } else {
            $min = $diff;
        }
    }
}

my $diff = $min + ($max - $min) / 2;
print "Approximate square root after $iterations iterations: $diff\n";

This of course is using floating point, but could easilly be addapted to fixed point. You can vary the accuracy by changing the iteration limit. Each iteration gets slightly more accurate than the one before.

eg: - find the square root of 9:

Approximate square root after 40 iterations: 2.99999999999955
   - or - 
Approximate square root after 10 iterations: 3.00048828125
   - or - 
Approximate square root after 5 iterations: 3.046875

If it had found the value 3 it would have stopped early of course.

Give it enough iterations and it should get it very accurate:

./sqrt.pl 0.00284
Square root found at 0.0532916503778969
59 iterations

It's the square root of the largest number you expect to deal with. — Majenko, Nov 04 '11 at 09:49

score 4 · Answer 2 · answered Nov 02 '11 at 10:34

Here are some ideas and routines from ascended transcendant master / Guru Scott Dattalo here .
That is of course a joke except the guru (Guru?) part. Scott is superb.

Relevant discussion. 2005 & PIC and some is C but may be of value.

Scott again - 2003

Two Masters !!!
Dattallo & Golovchenko.
A range of methods

score 3 · Answer 3 · answered Nov 02 '11 at 14:36

You didn't specify what you mean by "small value", or "approximation". So what I'm about to propose might not work, but here goes.

The easiest thing would be to make a look-up-table. Essentially a ROM where the address bus is the number that you want square-rooted, and the data output is the result. With a single BRAM, you could do a 9 bit in, 8 bit out LUT. Of course, more BRAM's will give you a bigger table.

(BRAM = The Xilinx term for a Block RAM, which can also be used as a ROM. Other FPGA's have similar things.)

If you want more precision than BRAM's will give you, you could do a simple linear interpolation of two LUT entries. For example, let's say that you want a 12-bit input, but you only have BRAM's for 10 bits. You take the top 10 bits of your input and look that up in the LUT. Add 1 to those 10 bits and look up that value too. You then do simple linear interpolation between the two results, using the bottom 2 bits to tell you the proportion of one value over the other. Of course this will only give you an approximation, but I think if you do the math you'll find that it just might be good enough.

This method is the least accurate with low-value numbers, but as the input goes to higher values the accuracy goes way up.

An optimization of the above method would be to use the BRAM's as a dual-port ROM. In this way you can read two values out without increasing the number of BRAM's used. This will also allow you to calculate a SQRT for every clock cycle, with some pipelining delays.

Incidentally, this method also works for SINE/COSINE too!

Small value means x approaching 0, that's why I'm interesting in the small value approximation of \sqrt{x}. — Ang Zhi Ping, Nov 04 '11 at 00:55
@angzhiping "Approaching zero" doesn't help. We need to know the range and accuracy. What you gave is half of the range and none of the accuracy. The end result is to know the number of input and output bits. Also important is the speed required: in terms of clock speed and clocks per sqrt. — , Nov 04 '11 at 02:48

score 3 · Answer 4 · answered Nov 03 '11 at 02:25

Try the following approach

If the number is negative, handle accordingly.
If the number is 0, return 0.
Otherwise:
normalize to a number in the range [1/4, 1]: count how many times k you have to multiply your number by 4 (x <<= 2 in C) until it is within the above range.
use an arbitrary approach (polynomial approximations, Newton's method for sqrt a[n] = (a[n-1]+k/a[n-1])/2, etc.) to calculate square root within this range
denormalize: shift right by k bits

score 0 · Answer 5 · answered Sep 25 '14 at 07:59

Try: improved guessing for 1st variable

Your number can be considered: A * 2^n
The 1st approximation is then: A * 2^(n/2)

Say you are using a 32-bit number, with 24 bits used to hold fractions. For numbers >1:
1. Count the number of bits used in the integer portion (N)
2. Halve this number (N' = N/2, ie. right shifted 1 bit)
3. Right shift the original number by N' : this is your 1st guess.

In this format, the smallest number you can have is 2^-24. The square root will ba about 2^-12. So, for numbers <1:
1. Count the number of "zero" bits in the fraction, until your reach a set bit (N)
2. Halve this number (N' = N/2, ie. right shifted 1 bit)
3. LEFT shift the original number by the revised count: this is your 1st guess.

Example:
0.0000 0000 0000 0000 1 [16 leading zeroes] approximates to: 0.0000 0000 1

Finally, if you still have trouble with small A: can you calculate 1/A?
If so, then invert your number, then try using the Inverse Square Root algorithm:
x' = 0.5x * (3 - Ax^2)

score 0 · Answer 6 · edited Nov 02 '12 at 19:03

0

Try \$x=(y+d)^2 \approx y^2+2dy\$ so let \$d=(x-y^2)/2y = (x/y-y) \gg 1\$ and next \$y=y+d. \$ If MSb is n from right, let first \$y=1 \ll (n/2)\$. Converges in <4 iterations.

edited Nov 02 '12 at 19:03

Nick Alexeev

38,181
17
100
237

answered Mar 14 '12 at 06:41

Byron

1

How to perform small value approximation for sqrt(x) on FPGA

6 Answers6