How is the square root function implemented? [closed]

Solution 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

Solution 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.

Solution 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.