Calculate an integer square root

From CodeCodex

Related content:

Contents

Algorithm

The primary algorithm used below to calculate the integer square root of a number n is a modified form of Newton's method for approximating roots (other correct algorithms may of course be used as well).

  1. Declare initial approximation
  2. Calculate next approximation
  3. Continue calculation until desired accuracy is achieved
  4. Compensate for error
    1. Error may occur when the square of the final approximation is greater than n. In given case, subtract one until the final approximation is less than n.

However, if you don't have a hardware multiply instruction, there is a faster basic algorithm.

Implementations

C

/*
Caution added by Martin L. Buchanan, mlb@backgroundtask.com, Wed 11/16/2005:

If number is the maximum unsigned int value, call it MAX_VAL, then the first 
evaluation of NEXT(n, number), with n == 1, produces an overflow when 
1 + MAX_VAL/1 is evaluated. For an unsigned type the overflow typically 
wraps around and yields zero as the macro result and zero as the 
overall function result.
*/
#include <stdlib.h>

#define NEXT(n, i)  (((n) + (i)/(n)) >> 1)

unsigned int isqrt(int number) {
  unsigned int n  = 1;
  unsigned int n1 = NEXT(n, number);

  while(abs(n1 - n) > 1) {
    n  = n1;
    n1 = NEXT(n, number);
  }
  while(n1*n1 > number)
    n1--;
  return n1;
}

But the following algorithm is many times faster, even when you have hardware multiplication and division.

unsigned long isqrt(unsigned long x)
{
    register unsigned long op, res, one;

    op = x;
    res = 0;

    /* "one" starts at the highest power of four <= than the argument. */
    one = 1 << 30;  /* second-to-top bit set */
    while (one > op) one >>= 2;

    while (one != 0) {
        if (op >= res + one) {
            op -= res + one;
            res += one << 1;  // <-- faster than 2 * one
        }
        res >>= 1;
        one >>= 2;
    }
    return res;
}

For an explanation as to how the previous and following algorithms work, see Integer Square Roots by Jack W. Crenshaw (http://www.embedded.com/98/9802fe2.htm) - specifically, Figure 2 of that article.
The fastest integer square root C algorithm yet is possibly below:

typedef unsigned char       uint8;
typedef unsigned short int  uint16;
typedef unsigned long int   uint32;

...

uint32  // OR uint16 OR uint8
isqrt32 (uint32 n) // OR isqrt16 ( uint16 n ) OR  isqrt8 ( uint8 n ) - respectively [ OR overloaded as isqrt (uint?? n) in C++ ]
{
    register uint32 // OR register uint16 OR register uint8 - respectively
        root, remainder, place;

    root = 0;
    remainder = n;
    place = 0x40000000; // OR place = 0x4000; OR place = 0x40; - respectively

    while (place > remainder)
        place = place >> 2;
    while (place)
    {
        if (remainder >= root + place)
        {
            remainder = remainder - root - place;
            root = root + (place << 1);
        }
        root = root >> 1;
        place = place >> 2;
    }
    return root;
}

However, this is completely processor-dependent and maybe compiler-dependent as well. Another very fast algorithm donated by Tristan Muntsinger (Tristan.Muntsinger@gmail.com) is below. It's best to test each one across an expected range of inputs to find the quickest one for your specific application.

unsigned int sqrt32(unsigned long n)
{
    unsigned int c = 0x8000;
    unsigned int g = 0x8000;

    for(;;) {
        if(g*g > n)
            g ^= c;
        c >>= 1;
        if(c == 0)
            return g;
        g |= c;
    }
}

C++

// Finds the integer square root of a positive number of any type
template <typename type>
type intSqrt (type remainder)
{
  if (remainder < 0) // if type is unsigned this will be ignored = no runtime
    return 0; // negative number ERROR

  type place = (type)1 << (sizeof (type) * 8 - 2); // calculated by precompiler = same runtime as: place = 0x40000000
  while (place > remainder)
    place /= 4; // optimized by complier as place >>= 2

  type root = 0;
  while (place)
  {
    if (remainder >= root+place)
    {
      remainder -= root+place;
      root += place * 2;
    }
    root /= 2;
    place /= 4;
  }
  return root;
}

The above does not work on signed integers larger than the initial value of "place", because "root += place * 2" will result in a negative value. That is, if the input is signed it must be less than 1/2 the maximum value of that type of integer.

To calculate the square root of the desired number,visit: http://bitsbyta.blogspot.com/2011/02/how-to-calculate-square-root-c.html For more related source codes visit: http://www.bitsbyta.blogspot.com/

Note that it is almost certainly faster to convert to floating point, use the built-in floating point sqrt, and then to convert back to int, than to use any of the above functions. But, beware loss of precision with large numbers. Using the GNU C libraries (libc6 v.2.12.1), the following is roughly 6.8x faster than the above function and returns the correct result for all numbers less than (1 << 24):

unsigned isqrt(unsigned value) {
  return static_cast<unsigned>(sqrt(static_cast<float>(value)));
}

Using the same libraries, integers up to (1 << 30) work if the double-precision sqrt is used instead, but this is only about 4x faster than the above algorithm. These cutoffs and timings are likely to vary depending on your CPU and OS, so be sure to test the output for all integers in the range that your program might use. But even with error-correcting code added, the built-in sqrt is significantly faster than any customized integer square root:

unsigned isqrt(unsigned value)
{
  unsigned result = static_cast<unsigned>(sqrt(static_cast<float>(value)));
  do { ++result; } while(result * result <= value);
  do { --result; } while(result * result > value);
  return result;
}

C#

// Finds the integer square root of a positive number
public static int Isqrt(int num) {
    if (0 == num) { return 0; }  // Avoid zero divide
    int n = (num / 2) + 1;       // Initial estimate, never low
    int n1 = (n + (num / n)) / 2;
    while (n1 < n) {
        n = n1;
        n1 = (n + (num / n)) / 2;
    } // end while
    return n;
} // end Isqrt()

Common Lisp

Common Lisp already includes isqrt. So, unless you want/have to, don't reinvent the wheel

* (isqrt 26)
5

Here's an implementation of the Newton's method (based on the C# implementation on this page), using tail-recursion.

 (defun integer-sqrt (num)
   (if (zerop num)
     0
     (labels ((next-approx (curr)
                           (truncate (/ (+ curr (truncate (/ num curr))) 2)))
              (int-sqrt (n n1)
                        (if (>= n1 n)
                          n
                          (int-sqrt n1 (next-approx n1)))))
       (let ((first-aprox (1+ (truncate (/ num 2)))))
         (int-sqrt first-aprox (next-approx first-aprox))))))

Erlang

isqrt(Num) when Num =< 0 -> 0;
isqrt(Num) ->
    isqrt(Num, 1, (1 + Num) div 2).

isqrt(Num, M, N) when abs(M - N) =< 1 ->
    if  N * N =< Num -> N;
        true         -> N - 1
    end;
isqrt(Num, _, N) ->
    isqrt(Num, N, (N + Num div N) div 2).

Haskell

 isqrt 0 = 0
 isqrt 1 = 1
 isqrt n = head $ dropWhile (\x -> x*x > n) $ iterate (\x -> (x + n `div` x) `div` 2) (n `div` 2)

Another solution:

 isqrt = floor . sqrt . fromIntegral

Java

The int version simply casts the result of StrictMath.sqrt to an int, giving us full hardware speed.

The long version uses a trick by Programmer Olathe to get exact results from StrictMath.sqrt even though doubles only have 52 bits precision: StrictMath.sqrt is guaranteed to give the exact result or the exact result plus one. This gets us very near full hardware speed.

longs have 64 bits and Java keeps only the most significant 52 bits when casting to double, so you can lose up to 12 bits.
If the input uses more than 52 bits, the distance between squares is at least (226)2 - (226 - 1)2 = 227 - 1. So, even if you lose the maximum of 12 bits, the distance between perfect squares is much greater than a 12-bit number, and so there cannot be more than one perfect square per "zone" of longs that cast to the same double.
This leaves two possibilities.
  1. There are no perfect squares in the zone, in which case the integer part of sqrt will correspond to the next lowest perfect square, which is great.
  2. There is only one perfect square in the zone, in which case the integer part of sqrt will correspond to that perfect square, and that's either great or, if our input is less than that perfect square, we simply need to subtract one.

The BigInteger version uses the long version to quickly get the square root of the most significant 64 bits, and then it fills in the remaining bits of the result from most to least significant, checking whether putting a 1 in that place would exceed the input and using 0 instead if that's the case. It is designed to minimize BigInteger object creation, creating less than two intermediate BigIntegers per output bit.

public class MathTools {
  // floorSqrt :: unsigned int -> unsigned int
  // Gives the exact floor of the square root of x, treated as unsigned.
  public static final int floorSqrt(final int x) {
    return (int) StrictMath.sqrt(x & 0xffffffffL);
  }
  
  // floorSqrt :: unsigned long -> unsigned int
  // Gives the exact floor of the square root of x, treated as unsigned.
  public static final int floorSqrt(final long x) {
    if ((x & 0xfff0000000000000L) == 0L) return (int) StrictMath.sqrt(x);
    final long result = (long) StrictMath.sqrt(2.0d*(x >>> 1));
    return result*result - x > 0L ? (int) result - 1 : (int) result;
  }
  
  // floorSqrt :: BigInteger -> BigInteger
  // Gives the exact floor of the square root of x, returning null (like Math.sqrt's NaN) if x is negative.
  public static final BigInteger floorSqrt(final BigInteger x) {
    if (x == null) return null;
    {  
      final int zeroCompare = x.compareTo(BigInteger.ZERO);  
      if (zeroCompare <  0) return null;  
      if (zeroCompare == 0) return BigInteger.ZERO;  
    }  
    
    int bit = Math.max(0, (x.bitLength() - 63) & 0xfffffffe); // last even numbered bit in first 64 bits
    BigInteger result = BigInteger.valueOf(floorSqrt(x.shiftRight(bit).longValue()) & 0xffffffffL);
    bit >>>= 1;
    result = result.shiftLeft(bit);
    while (bit != 0) {
      bit--;
      final BigInteger resultHigh = result.setBit(bit);
      if (resultHigh.multiply(resultHigh).compareTo(x) <= 0) result = resultHigh;
    }
    
    return result;
  }
}

Of course, if you just wanted a simple answer, here you go:

public class SquareRoot{
    public static void main(String args[]){
        int x=49;//store value in the variable to bee rooted
        double root=Math.sqrt(x); //take the root of x and store it
        System.out.println(root); //print the root
    }//end main method
}//end class SquareRoot

Also: you can typecast root into an int, meaning that you turn root into an int. like this:

public class SquareRoot{
    public static void main(String args[]){
        int x=49;
        int root= (int) Math.sqrt(x); //take the root of x , typecast it, then store it
        System.out.println(root); //print the root
    }//end main method
}//end class SquareRoot

OCaml

Recursive version:

let rec isqrt n =
  if n = 1 then 1
  else let n' = isqrt (n - 1) in
    (n' + (n / n')) / 2

Iterative version:

# let rec isqrt ?(xn=1) number =
   let xn1 = (xn + number / xn) lsr 1 in
   if xn1 = xn then xn1 else
     isqrt ~xn:xn1 number;;
val isqrt : ?xn:int -> int -> int = <fun>

For example:

# isqrt 12345;;
- : int = 111
# sqrt 12345.;;
- : float = 111.108055513540506

Perl

sub isqrt {
    my $num = shift;
    return 0 if 0 == $num;
    my $n = $num / 2 + 1;
    my $n1 = ($n + ($num / $n)) / 2;
    while ($n1 < $n) {
        $n = $n1;
        $n1 = ($n + ($num / $n)) / 2;
    };
    return $n;
};

Python

def isqrt(n):
    xn = 1
    xn1 = (xn + n/xn)/2
    while abs(xn1 - xn) > 1:
        xn = xn1
        xn1 = (xn + n/xn)/2
    while xn1*xn1 > n:
        xn1 -= 1
    return xn1

Python 3.1

This function uses int.bit_length() (available in Python 3.1) to get a better first approximation. It's much faster for large arguments. If bit_length isn't available, you can fake it by mucking around with bin().

def isqrt(n):
    'isqrt(n)\n\nReturn floor(sqrt(n)).'

    if not isinstance(n, int):
        raise TypeError('an int is required')
    if n < 0:
        raise ValueError('math domain error')

    guess = (n >> n.bit_length() // 2) + 1
    result = (guess + n // guess) // 2
    while abs(result - guess) > 1:
        guess = result
        result = (guess + n // guess) // 2
    while result * result > n:
        result -= 1
    return result

Ruby

def isqrt(square)
    # Check our input
    square = square.to_i # force integer
    return 0 if square == 0 # quick exit
    raise RangeError if square < 0
    # Actual computation
    n = iter(1, square)
    n1 = iter(n, square)
    n1, n = iter(n1, square), n1 until n1 >= n - 1
    n1 = n1 - 1 until n1*n1 <= square
    return n1
end
def iter(n, square) (n+(square/n))/2 end

Scheme

Recursive version:

(define (isqrt n)
  (if (= n 1)
      1
      (let ((xn (isqrt (- n 1))))
        (quotient (+ xn (quotient n xn)) 2))))

Iterative version:

(define (isqrt n)
  (let loop ((xn 1))
    (let ((xn1 (quotient (+ xn (quotient n xn)) 2)))
      (if (= xn1 xn)
          xn1
          (loop xn1)))))

Tcl

namespace import ::tcl::mathfunc::*
::tcl::mathfunc::isqrt 26

Seed7

const func integer: isqrt (in var integer: number) is func
  result
    var integer: result is 0;
  local
    var integer: res2 is 0;
  begin
    if number > 0 then
      res2 := number;
      repeat
        result := res2;
        res2 := (result + number div result) div 2;
      until result <= res2;
    end if;
  end func;

Zsh

Using mathfunc module

zmodload -i zsh/mathfunc
sqrt(n)

Newton's method implementation (based on the Perl implementation on this page)

zsqrt() {
	local -i num
	local -F n n1
	num=n
	((n = num / 2 + 1))
	((n1 = (n + (num / n)) / 2))
	while [ "${n1%.*}" -lt  "${n%.*}" ]
	do
		n=$n1
		((n1 = (n + (num / n)) / 2))
	done
	print $n1
}