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.
Answer from user448810 on Stack OverflowNote: 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.
Update: Python 3.8 has a math.isqrt function in the standard library!
I benchmarked every (correct) function here on both small (0…222) and large (250001) 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 largeint(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.)
Videos
» pip install isqrt
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) 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 are required on -25, otherwise it's parsed as -(25**.5) because exponentiation is more tightly binding than 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 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 let's try a big number:
>>> n = 10**30
>>> x = n**2
>>> root = x**.5
>>> root == n
False
>>> root - n # 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
SymPy
Depending on your goal, it might be a good idea to delay the calculation of square roots for as long as possible. SymPy might help.
SymPy is a Python library for symbolic mathematics.
import sympy
sympy.sqrt(2)
# => sqrt(2)
This doesn't seem very useful at first.
But sympy can give more information than floats or Decimals:
sympy.sqrt(8) / sympy.sqrt(27)
# => 2*sqrt(6)/9
Also, no precision is lost. (√2)² is still an integer:
s = sympy.sqrt(2)
s**2
# => 2
type(s**2)
#=> <class 'sympy.core.numbers.Integer'>
In comparison, floats and Decimals would return a number which is very close to 2 but not equal to 2:
(2**0.5)**2
# => 2.0000000000000004
from decimal import Decimal
(Decimal('2')**Decimal('0.5'))**Decimal('2')
# => Decimal('1.999999999999999999999999999')
Sympy also understands more complex examples like the Gaussian integral:
from sympy import Symbol, integrate, pi, sqrt, exp, oo
x = Symbol('x')
integrate(exp(-x**2), (x, -oo, oo))
# => sqrt(pi)
integrate(exp(-x**2), (x, -oo, oo)) == sqrt(pi)
# => True
Finally, if a decimal representation is desired, it's possible to ask for more digits than will ever be needed:
sympy.N(sympy.sqrt(2), 1_000_000)
# => 1.4142135623730950488016...........2044193016904841204