# Difference between revisions of "Calculate an integer square root"

### From CodeCodex

(→Ruby) |
m (→Java: +exact) |
||

(19 intermediate revisions by one user not shown) | |||

Line 166: | Line 166: | ||

=== Java === | === Java === | ||

− | + | A very fast algorithm by Programmer Olathe that takes an unsigned (yes, unsigned) int or long and returns an unsigned int. Uses the floating-point sqrt method to get a fairly precise result very quickly and then corrects it to an exact result. | |

− | + | ||

− | + | A double has 52-bit precision, so if the input uses 52 bits or less, use the floating-point result directly. | |

− | + | ||

− | + | ||

− | + | If the input uses more than 52 bits, the distance between squares is at least (2^26)^2 - (2^26 - 1)^2 = 2^27 - 1. So, even if the number uses all 64 bits, and 12 bits are lost when converting to double, the distance between perfect squares is much more than twice 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. | |

− | + | ||

− | + | # 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. | |

− | + | # 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. | |

− | + | ||

− | + | <pre class="java"> | |

− | + | public class MathTools { | |

− | return | + | // floorSqrt :: unsigned long -> 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); } | ||

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

} | } | ||

} | } |

## Revision as of 01:28, 3 December 2010

## Contents |

## Algorithm

The algorithm used in this source to calculate the integer square root of a number n is a modified form of Newton's method for approximating roots.

- Declare initial approximation
- Calculate next approximation
- Continue calculation until desired accuracy is achieved
- Compensate for error
- 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; }

### 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()

### Java

A very fast algorithm by Programmer Olathe that takes an unsigned (yes, unsigned) int or long and returns an unsigned int. Uses the floating-point sqrt method to get a fairly precise result very quickly and then corrects it to an exact result.

A double has 52-bit precision, so if the input uses 52 bits or less, use the floating-point result directly.

If the input uses more than 52 bits, the distance between squares is at least (2^26)^2 - (2^26 - 1)^2 = 2^27 - 1. So, even if the number uses all 64 bits, and 12 bits are lost when converting to double, the distance between perfect squares is much more than twice 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.

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

public class MathTools { // floorSqrt :: unsigned long -> 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); } 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; } }

### 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

### 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)))))

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

### 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))))))

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