Integer square root in python

Question or problem about Python programming:

Is there an integer square root somewhere in python, or in standard libraries? I want it to be exact (i.e. return an integer), and bark if there’s no solution.

At the moment I rolled my own naive one:

def isqrt(n):
    i = int(math.sqrt(n) + 0.5)
    if i**2 == n:
        return i
    raise ValueError('input was not a perfect square')

But it’s ugly and I don’t really trust it for large integers. I could iterate through the squares and give up if I’ve exceeded the value, but I assume it would be kinda slow to do something like that. Also I guess I’d probably be reinventing the wheel, something like this must surely exist in python already…

How to solve the problem:

Solution 1:

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.

Solution 2:

Update: Python 3.8 has a math.isqrt function in the standard library!

I benchmarked every (correct) function here on both small (0…2²²) and large (2⁵⁰⁰⁰¹) inputs. The clear winners in both cases are gmpy2.isqrt suggested by mathmandan in first place, followed by Python 3.8’s math.isqrt in second, followed by the ActiveState recipe linked by NPE in third. The ActiveState recipe has a bunch of divisions that can be replaced by shifts, which makes it a bit faster (but still behind the native functions):

def isqrt(n):
    if n > 0:
        x = 1 << (n.bit_length() + 1 >> 1)
        while True:
            y = (x + n // x) >> 1
            if y >= x:
                return x
            x = y
    elif n == 0:
        return 0
    else:
        raise ValueError("square root not defined for negative numbers")

Benchmark results:

• gmpy2.isqrt() (mathmandan): 0.08 µs small, 0.07 ms large
• int(gmpy2.isqrt())*: 0.3 µs small, 0.07 ms large
• Python 3.8 math.isqrt: 0.13 µs small, 0.9 ms large
• ActiveState (optimized as above): 0.6 µs small, 17.0 ms large
• ActiveState (NPE): 1.0 µs small, 17.3 ms large
• castlebravo long-hand: 4 µs small, 80 ms large
• mathmandan improved: 2.7 µs small, 120 ms large
• martineau (with this correction): 2.3 µs small, 140 ms large
• nibot: 8 µs small, 1000 ms large
• mathmandan: 1.8 µs small, 2200 ms large
• castlebravo Newton’s method: 1.5 µs small, 19000 ms large
• user448810: 1.4 µs small, 20000 ms large

(* Since gmpy2.isqrt returns a gmpy2.mpz object, which behaves mostly but not exactly like an int, you may need to convert it back to an int for some uses.)

Solution 3:

Sorry for the very late response; I just stumbled onto this page. In case anyone visits this page in the future, the python module gmpy2 is designed to work with very large inputs, and includes among other things an integer square root function.

Example:

>>> import gmpy2
>>> gmpy2.isqrt((10**100+1)**2)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001L)
>>> gmpy2.isqrt((10**100+1)**2 - 1)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000L)

Granted, everything will have the "mpz" tag, but mpz's are compatible with int's:

>>> gmpy2.mpz(3)*4
mpz(12)

>>> int(gmpy2.mpz(12))
12

See my other answer for a discussion of this method's performance relative to some other answers to this question.

Download: https://code.google.com/p/gmpy/

Solution 4:

Long-hand square root algorithm

It turns out that there is an algorithm for computing square roots that you can compute by hand, something like long-division. Each iteration of the algorithm produces exactly one digit of the resulting square root while consuming two digits of the number whose square root you seek. While the "long hand" version of the algorithm is specified in decimal, it works in any base, with binary being simplest to implement and perhaps the fastest to execute (depending on the underlying bignum representation).

Because this algorithm operates on numbers digit-by-digit, it produces exact results for arbitrarily large perfect squares, and for non-perfect-squares, can produce as many digits of precision (to the right of the decimal place) as desired.

There are two nice writeups on the "Dr. Math" site that explain the algorithm:

• Square Roots in Binary
• Longhand Square Roots

And here's an implementation in Python:

def exact_sqrt(x):
    """Calculate the square root of an arbitrarily large integer. 

    The result of exact_sqrt(x) is a tuple (a, r) such that a**2 + r = x, where
    a is the largest integer such that a**2 <= x, and r is the "remainder".  If
    x is a perfect square, then r will be zero.

    The algorithm used is the "long-hand square root" algorithm, as described at
    http://mathforum.org/library/drmath/view/52656.html

    Tobin Fricke 2014-04-23
    Max Planck Institute for Gravitational Physics
    Hannover, Germany
    """

    N = 0   # Problem so far
    a = 0   # Solution so far

    # We'll process the number two bits at a time, starting at the MSB
    L = x.bit_length()
    L += (L % 2)          # Round up to the next even number

    for i in xrange(L, -1, -1):

        # Get the next group of two bits
        n = (x >> (2*i)) & 0b11

        # Check whether we can reduce the remainder
        if ((N - a*a) << 2) + n >= (a<<2) + 1:
            b = 1
        else:
            b = 0

        a = (a << 1) | b   # Concatenate the next bit of the solution
        N = (N << 2) | n   # Concatenate the next bit of the problem

    return (a, N-a*a)

You could easily modify this function to conduct additional iterations to calculate the fractional part of the square root. I was most interested in computing roots of large perfect squares.

I'm not sure how this compares to the "integer Newton's method" algorithm. I suspect that Newton's method is faster, since it can in principle generate multiple bits of the solution in one iteration, while the "long hand" algorithm generates exactly one bit of the solution per iteration.

Source repo: https://gist.github.com/tobin/11233492

Solution 5:

Here's a very straightforward implementation:

def i_sqrt(n):
    i = n.bit_length() >> 1    # i = floor( (1 + floor(log_2(n))) / 2 )
    m = 1 << i    # m = 2^i
    #
    # Fact: (2^(i + 1))^2 > n, so m has at least as many bits 
    # as the floor of the square root of n.
    #
    # Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2)
    # >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED.
    #
    while m*m > n:
        m >>= 1
        i -= 1
    for k in xrange(i-1, -1, -1):
        x = m | (1 << k)
        if x*x <= n:
            m = x
    return m

This is just a binary search. Initialize the value m to be the largest power of 2 that does not exceed the square root, then check whether each smaller bit can be set while keeping the result no larger than the square root. (Check the bits one at a time, in descending order.)

For reasonably large values of n (say, around 10**6000, or around 20000 bits), this seems to be:

• Faster than the Newton's method implementation described by user448810.
• Much, much slower than the gmpy2 built-in method in my other answer.
• Comparable to, but somewhat slower than, the Longhand Square Root described by nibot.

All of these approaches succeed on inputs of this size, but on my machine, this function takes around 1.5 seconds, while @Nibot's takes about 0.9 seconds, @user448810's takes around 19 seconds, and the gmpy2 built-in method takes less than a millisecond(!). Example:

>>> import random
>>> import timeit
>>> import gmpy2
>>> r = random.getrandbits
>>> t = timeit.timeit
>>> t('i_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. # This function
1.5102493192883117
>>> t('exact_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. # Nibot
0.8952787937686366
>>> t('isqrt(r(20000))', 'from __main__ import *', number = 5)/5. # user448810
19.326695976676184
>>> t('gmpy2.isqrt(r(20000))', 'from __main__ import *', number = 5)/5. # gmpy2
0.0003599147067689046
>>> all(i_sqrt(n)==isqrt(n)==exact_sqrt(n)[0]==int(gmpy2.isqrt(n)) for n in (r(1500) for i in xrange(1500)))
True

This function can be generalized easily, though it's not quite as nice because I don't have quite as precise of an initial guess for m:

def i_root(num, root, report_exactness = True):
    i = num.bit_length() / root
    m = 1 << i
    while m ** root < num:
        m <<= 1
        i += 1
    while m ** root > num:
        m >>= 1
        i -= 1
    for k in xrange(i-1, -1, -1):
        x = m | (1 << k)
        if x ** root <= num:
            m = x
    if report_exactness:
        return m, m ** root == num
    return m

However, note that gmpy2 also has an i_root method.

In fact this method could be adapted and applied to any (nonnegative, increasing) function f to determine an "integer inverse of f". However, to choose an efficient initial value of m you'd still want to know something about f.

Edit: Thanks to @Greggo for pointing out that the i_sqrt function can be rewritten to avoid using any multiplications. This yields an impressive performance boost!

def improved_i_sqrt(n):
    assert n >= 0
    if n == 0:
        return 0
    i = n.bit_length() >> 1    # i = floor( (1 + floor(log_2(n))) / 2 )
    m = 1 << i    # m = 2^i
    #
    # Fact: (2^(i + 1))^2 > n, so m has at least as many bits
    # as the floor of the square root of n.
    #
    # Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2)
    # >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED.
    #
    while (m << i) > n: # (m<>= 1 i -= 1 d = n - (m << i) # d = n-m^2 for k in xrange(i-1, -1, -1): j = 1 << k new_diff = d - (((m<<1) | j) << k) # n-(m+2^k)^2 = n-m^2-2*m*2^k-2^(2k) if new_diff >= 0: d = new_diff m |= j return m 

Note that by construction, the kth bit of m << 1 is not set, so bitwise-or may be used to implement the addition of (m<<1) + (1<<k). Ultimately I have (2*m*(2**k) + 2**(2*k)) written as (((m<<1) | (1<<k)) << k), so it's three shifts and one bitwise-or (followed by a subtraction to get new_diff). Maybe there is still a more efficient way to get this? Regardless, it's far better than multiplying m*m! Compare with above:

>>> t('improved_i_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. 0.10908999762373242 >>> all(improved_i_sqrt(n) == i_sqrt(n) for n in xrange(10**6)) True 

Solution 6:

One option would be to use the decimal module, and do it in sufficiently-precise floats:

import decimal def isqrt(n): nd = decimal.Decimal(n) with decimal.localcontext() as ctx: ctx.prec = n.bit_length() i = int(nd.sqrt()) if i**2 != n: raise ValueError('input was not a perfect square') return i 

which I think should work:

>>> isqrt(1) 1 >>> isqrt(7**14) == 7**7 True >>> isqrt(11**1000) == 11**500 True >>> isqrt(11**1000+1) Traceback (most recent call last): File "", line 1, in isqrt(11**1000+1) File "", line 10, in isqrt raise ValueError('input was not a perfect square') ValueError: input was not a perfect square 

Solution 7:

Seems like you could check like this:

if int(math.sqrt(n))**2 == n: print n, 'is a perfect square' 

Update:

As you pointed out the above fails for large values of n. For those the following looks promising, which is an adaptation of the example C code, by Martin Guy @ UKC, June 1985, for the relatively simple looking binary numeral digit-by-digit calculation method mentioned in the Wikipedia article Methods of computing square roots:

from math import ceil, log def isqrt(n): res = 0 bit = 4**int(ceil(log(n, 4))) if n else 0 # smallest power of 4 >= the argument while bit: if n >= res + bit: n -= res + bit res = (res >> 1) + bit else: res >>= 1 bit >>= 2 return res if __name__ == '__main__': from math import sqrt # for comparison purposes for i in range(17)+[2**53, (10**100+1)**2]: is_perfect_sq = isqrt(i)**2 == i print '{:21,d}: math.sqrt={:12,.7G}, isqrt={:10,d} {}'.format( i, sqrt(i), isqrt(i), '(perfect square)' if is_perfect_sq else '') 

Output:

 0: math.sqrt= 0, isqrt= 0 (perfect square) 1: math.sqrt= 1, isqrt= 1 (perfect square) 2: math.sqrt= 1.414214, isqrt= 1 3: math.sqrt= 1.732051, isqrt= 1 4: math.sqrt= 2, isqrt= 2 (perfect square) 5: math.sqrt= 2.236068, isqrt= 2 6: math.sqrt= 2.44949, isqrt= 2 7: math.sqrt= 2.645751, isqrt= 2 8: math.sqrt= 2.828427, isqrt= 2 9: math.sqrt= 3, isqrt= 3 (perfect square) 10: math.sqrt= 3.162278, isqrt= 3 11: math.sqrt= 3.316625, isqrt= 3 12: math.sqrt= 3.464102, isqrt= 3 13: math.sqrt= 3.605551, isqrt= 3 14: math.sqrt= 3.741657, isqrt= 3 15: math.sqrt= 3.872983, isqrt= 3 16: math.sqrt= 4, isqrt= 4 (perfect square) 9,007,199,254,740,992: math.sqrt=9.490627E+07, isqrt=94,906,265 100,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,020,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001: math.sqrt= 1E+100, isqrt=10,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001 (perfect square) 

Solution 8:

Your function fails for large inputs:

In [26]: isqrt((10**100+1)**2) ValueError: input was not a perfect square 

There is a recipe on the ActiveState site which should hopefully be more reliable since it uses integer maths only. It is based on an earlier StackOverflow question: Writing your own square root function

Solution 9:

Floats cannot be precisely represented on computers. You can test for a desired proximity setting epsilon to a small value within the accuracy of python's floats.

def isqrt(n): epsilon = .00000000001 i = int(n**.5 + 0.5) if abs(i**2 - n) < epsilon: return i raise ValueError('input was not a perfect square') 

Solution 10:

I have compared the different methods given here with a loop:

for i in range (1000000): # 700 msec r=int(123456781234567**0.5+0.5) if r**2==123456781234567:rr=r else:rr=-1 

finding that this one is fastest and need no math-import. Very long might fail, but look at this

15241576832799734552675677489**0.5 = 123456781234567.0 

Hope this helps!