How is the square root function implemented?

Question

How is the square root function implemented?

Answer 1

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(100 * 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

Answer 2

On Intel hardware, it's often implemented on top of the hardware SQRT instruction. Some libraries just use the result of that straight off, some may put it through a couple of rounds of Newton optimisation to make it more accurate in the corner cases.

Answer 3

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.

Answer 4

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$

Answer 5

sqrt(); function Behind the scenes.

It always checks for the mid-points in a graph. Example: sqrt(16)=4; sqrt(4)=2;

Now if you give any input inside 16 or 4 like sqrt(10)==?

It finds the mid point of 2 and 4 i.e = x ,then again it finds the mid point of x and 4 (It excludes lower bound in this input). It repeats this step again and again until it gets the perfect answer i.e sqrt(10)==3.16227766017 .It lies b/w 2 and 4.All this in-built function are created using calculus,differentiation and Integration.

Answer 6

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

Answer 7

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);}

Answer 8

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));
}

Answer 9

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;
}

Answer 10

To calculate the square root (without using inbuilt math.sqrt function):

SquareRootFunction.java

public class SquareRootFunction {

    public double squareRoot(double value,int decimalPoints)
    {
        int firstPart=0;


        /*calculating the integer part*/
        while(square(firstPart)<value)
        {
            firstPart++;            
        }

        if(square(firstPart)==value)
            return firstPart;
        firstPart--;

        /*calculating the decimal values*/
        double precisionVal=0.1;
        double[] decimalValues=new double[decimalPoints];
        double secondPart=0;

        for(int i=0;i<decimalPoints;i++)
        {
            while(square(firstPart+secondPart+decimalValues[i])<value)
            {
                decimalValues[i]+=precisionVal;
            }

            if(square(firstPart+secondPart+decimalValues[i])==value)
            {
                return (firstPart+secondPart+decimalValues[i]);
            }

            decimalValues[i]-=precisionVal;
            secondPart+=decimalValues[i];
            precisionVal*=0.1;
        }

        return(firstPart+secondPart);

    }


    public double square(double val)
    {
        return val*val;
    }

}

MainApp.java

import java.util.Scanner;

public class MainApp {

public static void main(String[] args) {

    double number;
    double result;
    int decimalPoints;
    Scanner in = new Scanner(System.in);

    SquareRootFunction sqrt=new SquareRootFunction();   
    System.out.println("Enter the number\n");               
    number=in.nextFloat();  

    System.out.println("Enter the decimal points\n");           
    decimalPoints=in.nextInt();

    result=sqrt.squareRoot(number,decimalPoints);

    System.out.println("The square root value is "+ result);

    in.close();

    }

}

Answer 11

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/

Answer 12

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));
    }
}

Answer 13

long long int floorSqrt(long long int x) 
{
    long long r = 0;
    while((long)(1<<r)*(long)(1<<r) <= x){
        r++;
    }
    r--;
    long long b = r -1;
    long long ans = 1 << r;
    while(b >= 0){
        if(((long)(ans|1<<b)*(long)(ans|1<<b))<=x){
            ans |= (1<<b);
        }
        b--;
    }
    return ans;
}

Answer 14

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

Answer 15

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