Square Root Function

[shiftLynx]

Yesterday, I was wondering over how the square root function works in most standard libraries. I've learnt from one of my maths classes the Newton-Raphson method of determining roots of functions, so I sat down and figured out a nice formula. Maybe this will be useful to somebody.

Firstly, let b be the square that we're trying to find the root of (where b > 0). Let x be the positive root of that function. Hence:

x^2 = b
.: b - x^2 = 0

Let f(x) = b - x^2
Hence f'(x) = -2x

Using Newton-Raphson, x_n+1 = x_n - (b - x_n^2) / -2x_n
=> x_n+1 = x_n + (b - x_n^2) / 2x_n
=> x_n+1 = x_n + (1/2)((b - x_n^2) / x_n)
=> x_n+1 = x_n + (1/2)((b / x_n) - x_n)
=> x_n+1 = (1/2)(2x_n + (b / x_n) - x_n)
=> x_n+1 = (1/2)((b / x_n) + x_n)

Represented nicely:



I've found that starting with a first estimate as b/2 and doing about 5 iterations provides a fairly accurate result for smallish squares (less than 200ish).

Of course, the advantage of using Newton Raphson over any other iterative process is that it's fast!

Here is some code (which is actually slower than the SQR function in FreeBASIC):

[syntax="FreeBASIC"]
Function sqrt(b!) As Single
If b! = 0.0 Then
sqrt = 0
Exit Function
Elseif b! < 0.0 Then
sqrt = -1
Exit Function
End If

x! = b! / 2.0

prev! = -1.0
While x! <> prev!
prev! = x!
x! = 0.5 * ((b! / x!) + x!)
Wend

sqrt = x!
End Function
[/syntax]

Hopefully this will be useful and/or interesting to somebody.

[DrV]

Yesterday, I was wondering over how the square root function works in most standard libraries.

I assume by 'standard libraries' you mean language runtimes? In that case, most (for the IBM PC compatible, anyway ) probably just use the x87 (FPU) FSQRT instruction. I don't know how that's implemented on the silicon, but it's probably algorithmically something similar to your code - a repeated approximation to the desired accuracy.

[shiftLynx]

most (for the IBM PC compatible, anyway ) probably just use the x87 (FPU) FSQRT instruction.

I didn't realise there was an FPU instruction for this. Cool.