Determining Floating Point Square Root

Question

How do I determine the square root of a floating point number? Is the Newton-Raphson method a good way? I have no hardware square root either. I also have no hardware divide (but I have implemented floating point divide).

If possible, I would prefer to reduce the number of divides as much as possible since they are so expensive.

Also, what should be the initial guess to reduce the total number of iterations???

Thank you so much!

Solution 1

When you use Newton-Raphson to compute a square-root, you actually want to use the iteration to find the reciprocal square root (after which you can simply multiply by the input--with some care for rounding--to produce the square root).

More precisely: we use the function f(x) = x^-2 - n. Clearly, if f(x) = 0, then x = 1/sqrt(n). This gives rise to the newton iteration:

x_(i+1) = x_i - f(x_i)/f'(x_i)
        = x_i - (x_i^-2 - n)/(-2x_i^-3)
        = x_i + (x_i - nx_i^3)/2
        = x_i*(3/2 - 1/2 nx_i^2)

Note that (unlike the iteration for the square root), this iteration for the reciprocal square root involves no divisions, so it is generally much more efficient.

I mentioned in your question on divide that you should look at existing soft-float libraries, rather than re-inventing the wheel. That advice applies here as well. This function has already been implemented in existing soft-float libraries.

---

Edit: the questioner seems to still be confused, so let's work an example: sqrt(612). 612 is 1.1953125 x 2^9 (or b1.0011001 x 2^9, if you prefer binary). Pull out the even portion of the exponent (9) to write the input as f * 2^(2m), where m is an integer and f is in the range [1,4). Then we will have:

sqrt(n) = sqrt(f * 2^2m) = sqrt(f)*2^m

applying this reduction to our example gives f = 1.1953125 * 2 = 2.390625 (b10.011001) and m = 4. Now do a newton-raphson iteration to find x = 1/sqrt(f), using a starting guess of 0.5 (as I noted in a comment, this guess converges for all f, but you can do significantly better using a linear approximation as an initial guess):

x_0 = 0.5
x_1 = x_0*(3/2 - 1/2 * 2.390625 * x_0^2)
    = 0.6005859...
x_2 = x_1*(3/2 - 1/2 * 2.390625 * x_1^2)
    = 0.6419342...
x_3 = 0.6467077...
x_4 = 0.6467616...

So even with a (relatively bad) initial guess, we get rapid convergence to the true value of 1/sqrt(f) = 0.6467616600226026.

Now we simply assemble the final result:

sqrt(f) = x_n * f = 1.5461646...
sqrt(n) = sqrt(f) * 2^m = 24.738633...

And check: sqrt(612) = 24.738633...

Obviously, if you want correct rounding, careful analysis needed to ensure that you carry sufficient precision at each stage of the computation. This requires careful bookkeeping, but it isn't rocket science. You simply keep careful error bounds and propagate them through the algorithm.

If you want to correct rounding without explicitly checking a residual, you need to compute sqrt(f) to a precision of 2p + 2 bits (where p is precision of the source and destination type). However, you can also take the strategy of computing sqrt(f) to a little more than p bits, square that value, and adjust the trailing bit by one if necessary (which is often cheaper).

sqrt is nice in that it is a unary function, which makes exhaustive testing for single-precision feasible on commodity hardware.

You can find the OS X soft-float sqrtf function on opensource.apple.com, which uses the algorithm described above (I wrote it, as it happens). It is licensed under the APSL, which may or not be suitable for your needs.

Solution 2

Probably (still) the fastest implementation for finding the inverse square root and the 10 lines of code that I adore the most.

It's based on Newton Approximation, but with a few quirks. There's even a great story around this.

Solution 3

Easiest to implement (you can even implement this in a calculator):

def sqrt(x, TOL=0.000001):
    y=1.0
    while( abs(x/y -y) > TOL ):
        y= (y+x/y)/2.0
    return y

This is exactly equal to newton raphson:

y(new) = y - f(y)/f'(y)

f(y) = y^2-x and f'(y) = 2y

Substituting these values:

y(new) = y - (y^2-x)/2y = (y^2+x)/2y = (y+x/y)/2

If division is expensive you should consider: http://en.wikipedia.org/wiki/Shifting_nth-root_algorithm .

Shifting algorithms:

Let us assume you have two numbers a and b such that least significant digit (equal to 1) is larger than b and b has only one bit equal to (eg. a=1000 and b=10). Let s(b) = log_2(b) (which is just the location of bit valued 1 in b).

Assume we already know the value of a^2. Now (a+b)^2 = a^2 + 2ab + b^2. a^2 is already known, 2ab: shift a by s(b)+1, b^2: shift b by s(b).

Algorithm:

Initialize a such that a has only one bit equal to one and a^2<= n < (2*a)^2. 
Let q=s(a).    
b=a
sqra = a*a

For i = q-1 to -10 (or whatever significance you want):
    b=b/2
    sqrab = sqra + 2ab + b^2
    if sqrab > n:
        continue
    sqra = sqrab
    a=a+b

n=612
a=10000 (16)

sqra = 256

Iteration 1:
    b=01000 (8) 
    sqrab = (a+b)^2 = 24^2 = 576
    sqrab < n => a=a+b = 24

Iteration 2:
    b = 4
    sqrab = (a+b)^2 = 28^2 = 784
    sqrab > n => a=a

Iteration 3:
    b = 2
    sqrab = (a+b)^2 = 26^2 = 676
    sqrab > n => a=a

Iteration 4:
    b = 1
    sqrab = (a+b)^2 = 25^2 = 625
    sqrab > n => a=a

Iteration 5:
    b = 0.5
    sqrab = (a+b)^2 = 24.5^2 = 600.25
    sqrab < n => a=a+b = 24.5

Iteration 6:
    b = 0.25
    sqrab = (a+b)^2 = 24.75^2 = 612.5625
    sqrab < n => a=a


Iteration 7:
    b = 0.125
    sqrab = (a+b)^2 = 24.625^2 = 606.390625
    sqrab < n => a=a+b = 24.625

and so on.

Solution 4

A good approximation to square root on the range [1,4) is

def sqrt(x):
  y = x*-0.000267
  y = x*(0.004686+y)
  y = x*(-0.034810+y)
  y = x*(0.144780+y)
  y = x*(-0.387893+y)
  y = x*(0.958108+y)
  return y+0.315413

Normalise your floating point number so the mantissa is in the range [1,4), use the above algorithm on it, and then divide the exponent by 2. No floating point divisions anywhere.

With the same CPU time budget you can probably do much better, but that seems like a good starting point.

View more solutions

Comments

[amit]

this is basically a variation of the crossing method, and it is expected to converge much slower than newton raphson

[Veridian]

@ElKamina While I have implemented floating point divide, I would like to reduce the number of times I have to execute it since it is extremely expensive. The processor I am using has no hardware floating point divide unit.

[ElKamina]

@amit This method is equivalent to newton raphson method.

[ElKamina]

@starbox What is the efficiency of your multiplication compared to addition or subtraction? That is, how many additions is equal to multiplication? Does it depend on values that are being multiplied?

[ElKamina]

What is the advantage of using f(x) = x^-2 - n versus f(x) = x^2 - n ??

[Veridian]

@ElKamina Multiplication (in floating point) is slightly more efficient than my addition. I haven't counted it out yet. I'm going to try the answer above and get back to you if it doesn't pan out. Thanks

[Stephen Canon]

@ElKamina: As explained in my answer, using x^-2 - n avoids division in the iteration, making it much more efficient on typical hardware, while preserving quadratic convergence.

[Stephen Canon]

@starbox: I believe you can get the libgcc and compiler_rt sources from the public repositories. compiler_rt is also browsable on the web: llvm.org/viewvc/llvm-project/compiler-rt/trunk/lib/…

[Veridian]

@ElKamina, could you elaborate with an example of how to use the shifting n-th root algorithm for floating point. I read that page but it didn't look easy to implement. I have hardware shifters btw.

[ElKamina]

@starbox take a look at the pseudo code. Let me know if something is unclear.

[Veridian]

@StephenCanon, What do you mean a range reduction so that n is in [1,4). I have a 32-bit mantissa, I could just change the exponent to -1 and the number is guaranteed to be in the range of 0.5 to 1. Is that what you mean? But then the number doesn't converge. I tried using n = 612, and I need a starting guess around 0.01 for the algorithm to converge. If you could start with n = 612 as an example, that would work well for me. Thanks man.

[Veridian]

how do you normalize the mantissa so that it is in the range [1,4]? Please use x = 612 as an example.

[Veridian]

@ElKamina, how do i normalize to between [1,4), i don't understand that notation, can you provide an example with n = 612?

[Stephen Canon]

@starbox: 612 = 2.390625 x 2^8, so sqrt(612) = sqrt(2.390625) x sqrt(2^8) = sqrt(2.390625) x 2^4. Pull out the even portion of the exponent, and you're left with a value in the range [1,4).

[Stephen Canon]

@starbox: a linear approximation like 1.1 - 1/6 * f. (The 1/6 is optimal, and you can solve for the optimal first term; I can't remember the exact value off the top of my head, but it's very nearly 1.1)

[ElKamina]

@starbox I didn't ask you to normalize between [1,4). You are perhaps confusing my answer with another one. See my working solution for 612.

[sigfpe]

If you're using IEEE floats you can easily pick off the exponent. If it's even then the floating point number is if the form a*2^(2n) where a is in [1,2). So use the above algorithm to give the result sqrt(a)*2^n. If the exponent is odd we have a*2^(2n+1)=(2a)*2^(2n). This time a is in [2,4) and we use the above algorithm to give: sqrt(2*a)*2^n.

[Veridian]

@StephenCanon, just for clarification, how did you determine the first optimal term?

[Stephen Canon]

@starbox: it was a while ago, but I believe that I symbolically solved for the value that minimized the approximation error (some basic analysis will show this occurs when the maximum approximation error on the interior of the interval (which will occur where the derivative of f(x) = 1/x^2 is equal to -1/6) is equal to the approximation error at the endpoints of the interval. Or you can just plot the error and fuss with the constant until the error is symmetric.

[Veridian]

@StephenCanon, I apologize, but I still don't follow how you determined the initial guess. Is there a paper you could refer me to that proves this is the best initial guess? Thank you for your time.

[Stephen Canon]

@starbox: I don't think there is any such paper. Finding the linear approximation that minimizes the approximation error of the initial guess requires only high-school level mathematics (not what you're likely to find in the literature). To claim that it is "optimal" would be rather more subtle (indeed, it may not be), but it is easy to show that it is "good enough".

[Veridian]

@StephenCanon, since the square root will be on the range from 0.5 to 1, wouldn't 0.75 be optimal? But I guess you want change your initial guess using information from how large x^2 is. I know it requires high school level mathematics, but I am not following your steps. You stated the maximum approximation error is where the derivative of f(x) = x^-2 is equal to -1/6, but I really don't understand your intuition here.

[Stephen Canon]

@starbox: we're looking for the optimal linear approximation, not constant. It will have the form a(x) = ax + b, and will have the property that the approximation error at the endpoints will be equal to the maximum approximation error on the interior of the interval [1,4). Since the error at the endpoints is equal, the slope of the line must be the same as the slope of the line connecting (1,1) to (4,0.5). So we know that the approximation will have the form a(x) = -x/6 + b, and we need to solve for b.

[treehouse]

Link was invalid. Probably moved to beyond3d.com/content/articles/8 ?

[Adam Hyland]

It's actually a different blog post, earlier by 2 years. web.archive.org/web/20111126121744/https://codemaestro.com/… is an archive copy from 2011.