Integer Square Root in Python

Integer square root in python

Note: There is now math.isqrt in stdlib, available since Python 3.8.

Newton's method works perfectly well on integers:

def isqrt(n):
x = n
y = (x + 1) // 2
while y < x:
x = y
y = (x + n // x) // 2
return x

This returns the largest integer x for which x * x does not exceed n. If you want to check if the result is exactly the square root, simply perform the multiplication to check if n is a perfect square.

I discuss this algorithm, and three other algorithms for calculating square roots, at my blog.

How do I calculate square root in Python?

Option 1: math.sqrt()

The math module from the standard library has a sqrt function to calculate the square root of a number. It takes any type that can be converted to float (which includes int) as an argument and returns a float.

>>> import math
>>> math.sqrt(9)
3.0

Option 2: Fractional exponent

The power operator (**) or the built-in pow() function can also be used to calculate a square root. Mathematically speaking, the square root of a equals a to the power of 1/2.

The power operator requires numeric types and matches the conversion rules for binary arithmetic operators, so in this case it will return either a float or a complex number.

>>> 9 ** (1/2)
3.0
>>> 9 ** .5 # Same thing
3.0
>>> 2 ** .5
1.4142135623730951

(Note: in Python 2, 1/2 is truncated to 0, so you have to force floating point arithmetic with 1.0/2 or similar. See Why does Python give the "wrong" answer for square root?)

This method can be generalized to nth root, though fractions that can't be exactly represented as a float (like 1/3 or any denominator that's not a power of 2) may cause some inaccuracy:

>>> 8 ** (1/3)
2.0
>>> 125 ** (1/3)
4.999999999999999

Edge cases

Negative and complex

Exponentiation works with negative numbers and complex numbers, though the results have some slight inaccuracy:

>>> (-25) ** .5  # Should be 5j
(3.061616997868383e-16+5j)
>>> 8j ** .5 # Should be 2+2j
(2.0000000000000004+2j)

Note the parentheses on -25! Otherwise it's parsed as -(25**.5) because exponentiation is more tightly binding than unary negation.

Meanwhile, math is only built for floats, so for x<0, math.sqrt(x) will raise ValueError: math domain error and for complex x, it'll raise TypeError: can't convert complex to float. Instead, you can use cmath.sqrt(x), which is more more accurate than exponentiation (and will likely be faster too):

>>> import cmath
>>> cmath.sqrt(-25)
5j
>>> cmath.sqrt(8j)
(2+2j)

Precision

Both options involve an implicit conversion to float, so floating point precision is a factor. For example:

>>> n = 10**30
>>> x = n**2
>>> root = x**.5
>>> n == root
False
>>> n - root # how far off are they?
0.0
>>> int(root) - n # how far off is the float from the int?
19884624838656

Very large numbers might not even fit in a float and you'll get OverflowError: int too large to convert to float. See Python sqrt limit for very large numbers?

Other types

Let's look at Decimal for example:

Exponentiation fails unless the exponent is also Decimal:

>>> decimal.Decimal('9') ** .5
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: unsupported operand type(s) for ** or pow(): 'decimal.Decimal' and 'float'
>>> decimal.Decimal('9') ** decimal.Decimal('.5')
Decimal('3.000000000000000000000000000')

Meanwhile, math and cmath will silently convert their arguments to float and complex respectively, which could mean loss of precision.

decimal also has its own .sqrt(). See also calculating n-th roots using Python 3's decimal module

Calculating square root using only integer math in python

Note that the largest possible sum is 3*1024**2, so the largest possible square root is 1773 (floor - or 1774 rounded).

So you could simply take 0 as a starting guess, and repeatedly add 1 until the square exceeds the sum. That can't take more than about 1770 iterations.

Of course that's probably too slow. A straightforward binary search can cut that to 11 iterations, and doesn't require division (I'm assuming the MCU can shift right by 1 bit, which is the same as floor-division by 2).

EDIT

Here's some code, for a binary search returning the floor of the true square root:

def isqrt(n):
if n <= 1:
return n
lo = 0
hi = n >> 1
while lo <= hi:
mid = (lo + hi) >> 1
sq = mid * mid
if sq == n:
return mid
elif sq < n:
lo = mid + 1
result = mid
else:
hi = mid - 1
return result

To check, run:

from math import sqrt
assert all(isqrt(i) == int(sqrt(i)) for i in range(3*1024**2 + 1))

That checks all possible inputs given what you said - and since binary search is notoriously tricky to get right in all cases, it's good to check every case! It doesn't take long on a "real" machine ;-)

PROBABLY IMPORTANT

To guard against possible overflow, and speed it significantly, change the initialization of lo and hi to this:

    hi = 1
while hi * hi <= n:
hi <<= 1
lo = hi >> 1

Then the runtime becomes proportional to the number of bits in the result, greatly speeding smaller results. Indeed, for sloppy enough definitions of "close", you could stop right there.

FOR POSTERITY ;-)

Looks like the OP doesn't actually need square roots at all. But for someone who may, and can't afford division, here's a simplified version of the code, also removing multiplications from the initialization. Note: I'm not using .bit_length() because lots of deployed Python versions don't support that.

def isqrt(n):
if n <= 1:
return n
hi, hisq = 2, 4
while hisq <= n:
hi <<= 1
hisq <<= 2
lo = hi >> 1
while hi - lo > 1:
mid = (lo + hi) >> 1
if mid * mid <= n:
lo = mid
else:
hi = mid
assert lo + 1 == hi
assert lo**2 <= n < hi**2
return lo

from math import sqrt
assert all(isqrt(i) == int(sqrt(i)) for i in range(3*1024**2 + 1))

square root of a number greater than 10^2000 in Python 3

The usual square root methods convert the parameter to a float value before doing the calculation. As you saw, this does not work well with very large integers.

So use a function that is designed to work on arbitrarily large integers. Here is one, guaranteed to return correct integer part of the square root of any positive integer. This function drops the fractional part of the result, which may or may not be what you want. Since this function uses iteration it is also slower than the built-in square root routines. The Decimal module works on larger integers than the built-in routines but the precision of the values must be defined in advance so it does not work on arbitrarily large values.

import math

_1_50 = 1 << 50 # 2**50 == 1,125,899,906,842,624

def isqrt(x):
"""Return the integer part of the square root of x, even for very
large integer values."""
if x < 0:
raise ValueError('square root not defined for negative numbers')
if x < _1_50:
return int(math.sqrt(x)) # use math's sqrt() for small parameters
n = int(x)
if n <= 1:
return n # handle sqrt(0)==0, sqrt(1)==1
# Make a high initial estimate of the result (a little lower is slower!!!)
r = 1 << ((n.bit_length() + 1) >> 1)
while True:
newr = (r + n // r) >> 1 # next estimate by Newton-Raphson
if newr >= r:
return r
r = newr

Why does Python give the wrong answer for square root? What is integer division in Python 2?

In Python 2, sqrt=x**(1/2) does integer division. 1/2 == 0.

So x(1/2) equals x(0), which is 1.

It's not wrong, it's the right answer to a different question.

If you want to calculate the square root without an import of the math module, you'll need to use x**(1.0/2) or x**(1/2.). One of the integers needs to be a floating number.

Note: this is not the case in Python 3, where 1/2 would be 0.5 and 1//2 would instead be integer division.

Integer Square root of a number

The piece that is missing from the description is the so-called base case of the recursion. It is trivial, but necessary to specify: the integer square root of 0 is 0. By repeatedly recursing with a value that is one fourth (integer divison) of the current value, you'll eventually get to that base case.

I'm not fluent in SML, but I believe it should be something like this:

fun intSquareRoot 0 = 0
| intSquareRoot n =
let
val m = n div 4
val i = intSquareRoot m
val k = 2 * i + 1
in
if k * k <= n then k
else k - 1
end;

Exactness of integer square root in Python

Using python 2.7 on a machine with the C type long as a 64 bit integer and C type double implemented as a 64 bit IEEE floating point number yields

>>> import math
>>> x = (2<<53) + 1
>>> int(math.sqrt(x*x)) == x
False

I cheated and chose a number that 64 bit IEEE floating point can't represent exactly (but python's integer type can), (2<<53) + 1. Python 2.7 computes x*x as a python 2.7 long integer. (Note: This is distinct from a C long; python 2.7 can represent 2<<600 as an integer, but C cannot.)

Speaking of 2<<600,

>>> import math
>>> x = 2<<600
>>> int(math.sqrt(x*x)) == x
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
OverflowError: long int too large to convert to float

list of integer square roots python

import math
terms = [1,2,3,4,5,6,7,9,10,11,12,13,16,24,36]

integer_roots = []

for i in terms:
tmp = math.sqrt(i)
if (tmp.is_integer()):
integer_roots.append(int(tmp))

print(integer_roots)

Finding square roots in surd form

I don't see a better way than by factoring the given numbers, summing the multiplicities of the prime factors, extracting the even parts of the multiplicities and forming the square root of the products.

E.g.

√(84.375)=√(2²3.7.3.5³)=√(2²3²5³7)=2.3.5√(5.7)=30√35


Related Topics



Leave a reply



Submit