How is the square root function implemented?
15 Answers
Simple implementation using Binary Search with C++
double root(double n){
// Max and min are used to take into account numbers less than 1
double lo = min(1, n), hi = max(1, n), mid;
// Update the bounds to be off the target by a factor of 10
while(100 * lo * lo < n) lo *= 10;
while(0.01 * hi * hi > n) hi *= 0.1;
for(int i = 0 ; i < 100 ; i++){
mid = (lo+hi)/2;
if(mid*mid == n) return mid;
if(mid*mid > n) hi = mid;
else lo = mid;
}
return mid;
}
Note that while loop is most common with the binary search but Personally I prefer using for when dealing with decimal numbers, it saves some special cases handling and gets pretty accurate result from small loops like that 1000 or even 500 (Both will give the same result for almost all numbers but just to be safe).
Edit: Check out this Wikipedia article for various -special purpose- methods specialized in computing the square root.
Edit 2: Apply updates suggested by @jorgbrown to fix the function in case of input less than 1. Also, apply optimization to make the bounds off the target root by a factor of 10
The FDLIBM (Freely Distributable LIBM) has a quite nice documented version of sqrt. e_sqrt.c.
The have one version which uses integer arithmetic and a recurrence formula modifying one bit at a time.
Another method uses Newton's method. It starts with some black magic and a lookup table to get the first 8 bits and then applies the recurrence formula
y_{i+1} = 1/2 * ( y_i + x / y_i)
where x is the number we started with. This is the Babylonian method of Heron's method. It dates back to Hero of Alexandra in the first centuary AD.
There is another method called the Fast inverse square root or reciproot. which uses some "evil floating point bit level hacking" to find the value of 1/sqrt(x). i = 0x5f3759df - ( i >> 1 ); It exploits the binary representation of a float using the mantisse and exponent. If our number x is (1+m) * 2^e, where m is the mantissa and e the exponent and the result y = 1/sqrt(x) = (1+n) * 2^f. Taking logs
lg(y) = - 1/2 lg(x)
f + lg(1+n) = -1/2 e - 1/2 lg(1+m)
So we see the exponent part of the result is -1/2 the exponent of the number. The black magic basically does a bitwise shift on the exponent and uses a linear approximation on the mantissa.
Once you have a good first approximation you can use Newton's methods to get a better result and finally some bit level work to fix the last digit.
This is an implementation of Newton's algorithm, see https://tour.golang.org/flowcontrol/8.
func Sqrt(x float64) float64 {
// let initial guess to be 1
z := 1.0
for i := 1; i <= 10; i++ {
z -= (z*z - x) / (2*z) // MAGIC LINE!!
fmt.Println(z)
}
return z
}
The following is a mathematical explanation of the magic line. Suppose you want to find the root of the polynomial $f(x) = x^2 - a$. By Newton's method, you could start with an initial guess $x_0 = 1$. The next guess is $x_1 = x_0 - f(x_0)/f'(x_0)$, where $f'(x) = 2x$. Therefore, your new guess is
$x_1 = x_0 - (x_0^2 - a)/2x_0$
Solutions so far have been mainly floating-point... and have also assumed that a divide instruction is available and fast.
Here's a simple straightforward routine that doesn't use FP or divide. Every line computes another bit in the result, except for the first if statement which speeds up the routine when the input is small.
constexpr unsigned int root(unsigned int x) {
unsigned int i = 0;
if (x >= 65536) {
if ((i + 32768) * (i + 32768) <= x) i += 32768;
if ((i + 16384) * (i + 16384) <= x) i += 16384;
if ((i + 8192) * (i + 8192) <= x) i += 8192;
if ((i + 4096) * (i + 4096) <= x) i += 4096;
if ((i + 2048) * (i + 2048) <= x) i += 2048;
if ((i + 1024) * (i + 1024) <= x) i += 1024;
if ((i + 512) * (i + 512) <= x) i += 512;
if ((i + 256) * (i + 256) <= x) i += 256;
}
if ((i + 128) * (i + 128) <= x) i += 128;
if ((i + 64) * (i + 64) <= x) i += 64;
if ((i + 32) * (i + 32) <= x) i += 32;
if ((i + 16) * (i + 16) <= x) i += 16;
if ((i + 8) * (i + 8) <= x) i += 8;
if ((i + 4) * (i + 4) <= x) i += 4;
if ((i + 2) * (i + 2) <= x) i += 2;
if ((i + 1) * (i + 1) <= x) i += 1;
return i;
}
Implementation in Python: The floor of the root value is the output of this function. Example: The square root of 8 is 2.82842..., this function will give output '2'
def mySqrt(x):
# return int(math.sqrt(x))
if x==0 or x==1:
return x
else:
start = 0
end = x
while (start <= end):
mid = int((start + end) / 2)
if (mid*mid == x):
return mid
elif (mid*mid < x):
start = mid + 1
ans = mid
else:
end = mid - 1
return ans
I'm making a sqrt function too, 100000000 iterations takes 14 seconds, still nothing compared to 1 seconds by sqrt
double mysqrt(double n)
{
double x = n;
int it = 4;
if (n >= 90)
{
it = 6;
}
if (n >= 5000)
{
it = 8;
}
if (n >= 20000)
{
it = 10;
}
if (n >= 90000)
{
it = 11;
}
if (n >= 200000)
{
it = 12;
}
if (n >= 900000)
{
it = 13;
}
if (n >= 3000000)
{
it = 14;
}
if (n >= 10000000)
{
it = 15;
}
if (n >= 30000000)
{
it = 16;
}
if (n >= 100000000)
{
it = 17;
}
if (n >= 300000000)
{
it = 18;
}
if (n >= 1000000000)
{
it = 19;
}
for (int i = 0; i < it; i++)
{
x = 0.5*(x+n/x);
}
return x;
}
But the fastest implementation is:
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
float mysqrt(float n) {return 1/Q_rsqrt(n);}
there is some thing called Babylonian method.
static float squareRoot(float n)
{
/*We are using n itself as
initial approximation This
can definitely be improved */
float x = n;
float y = 1;
// e decides the accuracy level
double e = 0.000001;
while(x - y > e)
{
x = (x + y)/2;
y = n/x;
}
return x;
}
for more information link: https://www.geeksforgeeks.org/square-root-of-a-perfect-square/
So, just in case there are no specifications about whether not to use the in-built ceil or round function, here is a recursive approach in Java to finding the square root of an unsigned number using the Newton-Raphson method.
public class FindSquareRoot {
private static double newtonRaphson(double N, double X, double oldX) {
if(N <= 0) return 0;
if (Math.round(X) == Math.ceil(oldX))
return X;
return newtonRaphson(N, X - ((X * X) - N)/(2 * X), X);
}
//Driver method
public static void main (String[] args) {
System.out.println("Square root of 48.8: " + newtonRaphson(48.8, 10, 0));
}
}
Formula: root(number, <root>, <depth>) == number^(root^(-depth))
Usage: root(number,<root>,<depth>)
Example: root(16,2) == sqrt(16) == 4
Example: root(16,2,2) == sqrt(sqrt(16)) == 2
Example: root(64,3) == 4
Implementation in C#:
static double root(double number, double root = 2f, double depth = 1f)
{
return Math.Pow(number, Math.Pow(root, -depth));
}
Following my solution in Golang.
package main
import (
"fmt"
)
func Sqrt(x float64) float64 {
z := 1.0 // initial guess to be 1
i := 0
for int(z*z) != int(x) { // until find the first approximation
// Newton root algorithm
z -= (z*z - x) / (2 * z)
i++
}
return z
}
func main() {
fmt.Println(Sqrt(8900009870))
}
Following a classic/common solution.
package main
import (
"fmt"
"math"
)
func Sqrt(num float64) float64 {
const DIFF = 0.0001 // To fix the precision
z := 1.0
for {
z1 := z - (((z * z) - num) / (2 * z))
// Return a result when the diff between the last execution
// and the current one is lass than the precision constant
if (math.Abs(z1 - z) < DIFF) {
break
}
z = z1
}
return z
}
func main() {
fmt.Println(Sqrt(94339))
}
For further information check here
Usage: root(number,root,depth)
Example: root(16,2) == sqrt(16) == 4
Example: root(16,2,2) == sqrt(sqrt(16)) == 2
Example: root(64,3) == 4
Implementation in C#:
double root(double number, double root, double depth = 1f)
{
return number ^ (root ^ (-depth));
}
Usage: Sqrt(Number,depth)
Example: Sqrt(16) == 4
Example: Sqrt(8,2) == sqrt(sqrt(8))
double Sqrt(double number, double depth = 1) return root(number,2,depth);
By: Imk0tter
Related