Calculate an integer square root
From CodeCodex
| Image:Lightbulb.png Related content: |
| <DPL>
category=Math shownamespace=false </DPL> |
Contents |
[edit] 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.
[edit] Implementations
[edit] 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?...
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;
}
[edit] 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;
}
[edit] 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()
[edit] Java
public class MathTools {
private static int next(int n, int i) {
return (n + i/n) >> 1;
}
public static int isqrt(int number) {
int n = 1;
int n1 = next(n, number);
while(Math.abs(n1 - n) > 1) {
n = n1;
n1 = next(n, number);
}
while((n1*n1) > number) {
n1 -= 1;
}
return n1;
}
}
[edit] Haskell
isqrt 1 = 1 isqrt n = (n' + n `div` n') `div` 2 where n' = isqrt (n-1)
Beware: this function is incorrect!
isqrt 1940449 *** Exception: stack overflow
isqrt n = head $ dropWhile (\x -> x*x > n) $ iterate (\x -> (x + n `div` x) `div` 2) (n `div` 2)
mucho betteros, but very ugly :P
and the best
isqrt = floor . sqrt . fromIntegral
[edit] 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
[edit] 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;
};
[edit] 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
[edit] 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
[edit] 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
[edit] 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)))))
[edit] 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;
[edit] 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))))))
[edit] 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
}
Categories: Assembly | C | C sharp | Haskell | Objective Caml | Java | Perl | Python | Seed7 | Ruby | Scheme | Common Lisp | Zsh | Math
