How do I calculate square root in Python?
Asked Answered
A

10

53

I need to calculate the square root of some numbers, for example √9 = 3 and √2 = 1.4142. How can I do it in Python?

The inputs will probably be all positive integers, and relatively small (say less than a billion), but just in case they're not, is there anything that might break?


Note: This is an attempt at a canonical question after a discussion on Meta about an existing question with the same title.

Related

Apulia answered 20/1, 2022 at 21:16 Comment(1)
Comments are not for extended discussion; this conversation has been moved to chat.Haricot
A
87

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

Apulia answered 20/1, 2022 at 21:16 Comment(5)
„Exponentiation works with negative numbers and complex numbers, though the results are very slightly off and I'm not sure why:“ That’s due to both being complex number operations/results and complex numbers not being a number line (from -inf to +inf) but rather a 2D plane (also -inf j and +inf j). Compare to how √x=1 has the solutions +1 and -1 - i.e. both „directions“ of the number line. Since complex numbers represent a plane, sqrt results are a circle on the complex plane. Picking one result on this circle is not numerically stable, hence some algorithms produce inaccurate results.Cutwork
@Mister Huh, cool! That begs the question though, why does cmath use a different algo?Apulia
@Mister Oh wait, if you square the results again, you don't get the input number back. Does that mean the actual problem is some sort of drop of precision?Apulia
@wjandrea: "why does cmath use a different algo?" <- Different from what? Are you asking why cmath.sqrt(z) doesn't just use z ** 0.5? If so, the answer is that general complex power is a more complicated algorithm (take complex log, scale, then take complex exp of the result) than square root, with more opportunities for loss of accuracy, and so cmath.sqrt(z) is likely to be both faster and more accurate than z ** 0.5. My recommendation would be to always use an explicit sqrt function or method rather than a general powering operation.Jorgan
@Mark Thanks, that makes a lot of sense. I tried cmath.exp(cmath.log(x)/2) and got results that were comparable in accuracy to x ** .5. I updated the answer accordingly. If there's anything you'd like to add, by all means. I was thinking maybe it could use a summary at the bottom saying, like, "x ** .5 works in general," [which, btw, is why I put it first] "but for more specific cases, use the .sqrt() function of the library you're using for optimal results and performance."Apulia
G
25

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
Gotthelf answered 21/1, 2022 at 14:38 Comment(0)
A
15

NumPy

>>> import numpy as np
>>> np.sqrt(25)
5.0
>>> np.sqrt([2, 3, 4])
array([1.41421356, 1.73205081, 2.        ])

docs

Negative

For negative reals, it'll return nan, so np.emath.sqrt() is available for that case.

>>> a = np.array([4, -1, np.inf])
>>> np.sqrt(a)
<stdin>:1: RuntimeWarning: invalid value encountered in sqrt
array([ 2., nan, inf])
>>> np.emath.sqrt(a)
array([ 2.+0.j,  0.+1.j, inf+0.j])

Another option, of course, is to convert to complex first:

>>> a = a.astype(complex)
>>> np.sqrt(a)
array([ 2.+0.j,  0.+1.j, inf+0.j])
Apulia answered 20/1, 2022 at 21:16 Comment(0)
J
5

Python's fractions module and its class, Fraction, implement arithmetic with rational numbers. The Fraction class doesn't implement a square root operation, because most square roots are irrational numbers. However, it can be used to approximate a square root with arbitrary accuracy, because a Fraction's numerator and denominator are arbitrary-precision integers.

The following method takes a positive number x and a number of iterations, and returns upper and lower bounds for the square root of x.

from fractions import Fraction

def sqrt(x, n):
    x = x if isinstance(x, Fraction) else Fraction(x)
    upper = x + 1
    for i in range(0, n):
        upper = (upper + x/upper) / 2
    lower = x / upper
    if lower > upper:
        raise ValueError("Sanity check failed")
    return (lower, upper)

See the reference below for details on this operation's implementation. It also shows how to implement other operations with upper and lower bounds (although there is apparently at least one error with the log operation there).

  • Daumas, M., Lester, D., Muñoz, C., "Verified Real Number Calculations: A Library for Interval Arithmetic", arXiv:0708.3721 [cs.MS], 2007.

Alternatively, using Python's math.isqrt, we can calculate a square root to arbitrary precision:

  • Square root of i within 1/2n of the correct value, where i is an integer:Fraction(math.isqrt(i * 2**(n*2)), 2**n).
  • Square root of i within 1/10n of the correct value, where i is an integer:Fraction(math.isqrt(i * 10**(n*2)), 10**n).
  • Square root of x within 1/2n of the correct value, where x is a multiple of 1/2n:Fraction(math.isqrt(x * 2**(n)), 2**n).
  • Square root of x within 1/10n of the correct value, where x is a multiple of 1/10n:Fraction(math.isqrt(x * 10**(n)), 10**n).

In the foregoing, i or x must be 0 or greater.

June answered 21/1, 2022 at 19:45 Comment(1)
This is similar to my method, except it terminates another way. I really like this one. Also, the fractions can be swapped for other similar types, i.e. fixed precision or others.Dolliedolloff
A
4

Newton's method

Most simple and accurate way to compute square root is Newton's method.

You have a number which you want to compute its square root (num) and you have a guess of its square root (estimate). Estimate can be any number bigger than 0, but a number that makes sense shortens the recursive call depth significantly.

new_estimate = (estimate + num/estimate) / 2

This line computes a more accurate estimate with those 2 parameters. You can pass new_estimate value to the function and compute another new_estimate which is more accurate than the previous one or you can make a recursive function definition like this.

def newtons_method(num, estimate):
    # Computing a new_estimate
    new_estimate = (estimate + num/estimate) / 2
    print(new_estimate)
    # Base Case: Comparing our estimate with built-in functions value
    if new_estimate == math.sqrt(num):
        return True
    else:
        return newtons_method(num, new_estimate)

For example we need to find 30's square root. We know that the result is between 5 and 6.

newtons_method(30,5)

number is 30 and estimate is 5. The result from each recursive calls are:

5.5
5.477272727272727
5.4772255752546215
5.477225575051661

The last result is the most accurate computation of the square root of number. It is the same value as the built-in function math.sqrt().


This answer was originally posted by gunesevitan, but is now deleted.

Apulia answered 20/1, 2022 at 21:16 Comment(6)
I'm always amazed by how simple this method is, and how fast it converges. It was known loooooong before Newton, BTW.Gotthelf
@EricDuminil Still, a Python implementation of Newton's method for sqrt will probably always be much slower than the builtin sqrt function and therefore not recommended for practical purposes. The answer would be improved if it would discuss when this implementation should be used (higher accuracy maybe even or faster for arbitrary roots or showcasing what math.sqrt likely does)?Wormseed
@Trilarion math.sqrt probably directly calls a processor's square root instruction. The only good reason I can think of for using something else is when you have unique requirements, e.g. looking for an integer square root.Marchall
You can swap types in this version too, and using fractions or fixed point numbers can fix the problem of this being a slower math.sqrt too.Dolliedolloff
I don't think it is a good idea to implement such a critical function in Python itself. I'd rather make this a C function and call that from Python.Newbill
@Newbill I don't think anyone's going to be using this solution for critical applications. I see it as a helpful learning tool at best.Apulia
D
4

Binary search

Disclaimer: this is for a more specialised use-case. This method might not be practical in all circumstances.

Benefits:

  • can find integer values (i.e. which integer is the root?)
  • no need to convert to float, so better precision (can be done that well too)

I personally implemented this one for a crypto CTF challenge (RSA cube root attack),where I needed a precise integer value.

The general idea can be extended to any other root.

def int_squareroot(d: int) -> tuple[int, bool]:
    """Try calculating integer squareroot and return if it's exact"""
    left, right = 1, (d+1)//2
    while left<right-1:
        x = (left+right)//2
        if x**2 > d:
            left, right = left, x
        else:
            left, right = x, right
    return left, left**2==d

EDIT:

As @wjandrea have also pointed out, **this example code can NOT compute **. This is a side-effect of the fact that it does not convert anything into floats, so no precision is lost. If the root is an integer, you get that back. If it's not, you get the biggest number whose square is smaller than your number. I updated the code so that it also returns a bool indicating if the value is correct or not, and also fixed an issue causing it to loop infinitely (also pointed out by @wjandrea). This implementation of the general method still works kindof weird for smaller numbers, but above 10 I had no problems with.

Overcoming the issues and limits of this method/implementation:

For smaller numbers, you can just use all the other methods from other answers. They generally use floats, which might be a loss of precision, but for small integers that should mean no problem at all. All of those methods that use floats have the same (or nearly the same) limit from this.

If you still want to use this method and get float results, it should be trivial to convert this to use floats too. Note that that will reintroduce precision loss, this method's unique benefit over the others, and in that case you can also just use any of the other answers. I think the newton's method version converges a bit faster, but I'm not sure.

For larger numbers, where loss of precision with floats come into play, this method can give results closer to the actual answer (depending on how big is the input). If you want to work with non-integers in this range, you can use other types, for example fixed precision numbers in this method too.

Edit 2, on other answers:

Currently, and afaik, the only other answer that has similar or better precision for large numbers than this implementation is the one that suggest SymPy, by Eric Duminil. That version is also easier to use, and work for any kind of number, the only downside is that it requires SymPy. My implementation is free from any huge dependencies if that is what you are looking for.

Dolliedolloff answered 21/1, 2022 at 17:12 Comment(10)
If d is not a perfect square, the result will always be lower than its actual root, right? There's an existing question about that: Integer square root in python. You might want to post this solution there as well. There's also How to find integer nth roots?.Apulia
@wjandrea: The above function seems to answer the question "what's the largest integer whose square fits in d?" It's similar to n // p and "how many times can p fit in d"?Gotthelf
Well done, BTW. It even works fine for int_squareroot(int('123456789'*100)**2).Gotthelf
@Apulia yes, this implementation does that. The general idea can be used with floats as well to achieve arbitrary precision (bound by float's precision), but then the other answer's methods are better in most terms.Dolliedolloff
This doesn't work for d=2: it never finishes. Even if it did finish, it'd return 1 (IIUC), which is not what the question is looking for: "√2 = 1.4142". I hate to downvote a potentially useful answer, but this just doesn't answer the question. You might want to post this on Integer square root in python instead, after fixing the algo.Apulia
It also doesn't work for d=1 or d=0.Apulia
After the edit, it doesn't work for d=4 -> (1, False). Why would you use this implementation instead of math.isqrt(), which always returns the correct result and is implemented in C? (for CPython). Either way, it's still not what the question is looking for, which is a real number result.Apulia
It also still doesn't work for d=0 -> (1, False). Maybe zero could be considered a special case, IDK, but math.isqrt() returns 0, which fits your criterion (left**2==d).Apulia
Yes, Newton's method is generally faster. Bisection increases the precision of the result by 1 bit per loop. Newton's method, starting with a reasonable first approximation, (roughly) doubles the precision per loop. With a bad 1st approximation, Newton's method degrades to bisection until the value gets close enough to the root. Python floats have 53 bits of precision.Misplace
@UncleDino my answer has similar or better precision (i.e., arbitrary precision)... while being correctly rounded (rather than truncated à la isqrt): https://mcmap.net/q/175991/-how-do-i-calculate-square-root-in-pythonOveruse
M
4

Arbitrary precision square root

This variation uses string manipulations to convert a string which represents a decimal floating-point number to an int, calls math.isqrt to do the actual square root extraction, and then formats the result as a decimal string. math.isqrt rounds down, so all produced digits are correct.

The input string, num, must use plain float format: 'e' notation is not supported. The num string can be a plain integer, and leading zeroes are ignored.

The digits argument specifies the number of decimal places in the result string, i.e., the number of digits after the decimal point.

from math import isqrt

def str_sqrt(num, digits):
    """ Arbitrary precision square root

        num arg must be a string
        Return a string with `digits` after
        the decimal point

        Written by PM 2Ring 2022.01.26
    """

    int_part , _, frac_part = num.partition('.')
    num = int_part + frac_part

    # Determine the required precision
    width = 2 * digits - len(frac_part)

    # Truncate or pad with zeroes
    num = num[:width] if width < 0 else num + '0' * width
    s = str(isqrt(int(num)))

    if digits:
        # Pad, if necessary
        s = '0' * (1 + digits - len(s)) + s
        s = f"{s[:-digits]}.{s[-digits:]}"
    return s

Test

print(str_sqrt("2.0", 30))

Output

1.414213562373095048801688724209

For small numbers of digits, it's faster to use decimal.Decimal.sqrt. Around 32 digits or so, str_sqrt is roughly the same speed as Decimal.sqrt. But at 128 digits, str_sqrt is 2.2× faster than Decimal.sqrt, at 512 digits, it's 4.3× faster, at 8192 digits, it's 7.4× faster.

Here's a live version running on the SageMathCell server.

Misplace answered 25/1, 2022 at 14:39 Comment(1)
Also see https://mcmap.net/q/46855/-integer-square-root-in-pythonMisplace
O
2

A correctly rounded, arbitrary-precision square root of any int, float, or Fraction, without looping or converting to Decimal

(Copied mostly from my original answer here)

Note: requires python 3.8+ as I'm using math.isqrt under the hood to work the magic.

def sqrt(x: Union[int, float, Fraction], precision: int = 53) -> Fraction:
    a, b = x.as_integer_ratio()
    d = a.bit_length() - b.bit_length()
    s = max(precision - 1 - (d-(b<<(d>0 and d)>a<<(d<0 and -d))>>1), 0)
    a <<= s << 1
    n0 = math.isqrt(a // b)
    n1 = n0 + 1
    return Fraction(n1 if n0 * n1 * b < a else n0, 1 << s)

More precisely, the following are guaranteed to hold:

  • |√(x) - sqrt(x, precision=p)| < 0.5ulpₚ(√(x))
  • sqrt(float(x), precision=53) == math.sqrt(x) as long as float(x) <= 2**106 (after which point, math.sqrt will be less precise).
  • sqrt(x * x) == x if x is an integer, avoiding this problem
  • More generally, sqrt(x, precision=p) will be correctly rounded to precision max(p, ⌊√(x)⌋.bit_length())
  • For the mathematically curious: d-(b<<(d>0 and d)>a<<(d<0 and -d))>>1 is equal to ⌊log₂(√(a/b))⌋ for a > 0.

If you need the result as a float instead of a Fraction, just do float(sqrt(x)) (although this may lose precision or overflow if the final result is too big for a float).

Note: if you are only operating on integers, there is a slightly simpler function, which is equivalent to the above for any integer or Fraction with a denominator of 1:

def sqrt_of_int(x: int, precision: int = 53) -> Fraction:
    s = max(precision - (x.bit_length() + 1 >> 1), 0)
    x <<= s << 1
    n0 = math.isqrt(x)
    n1 = n0 + 1
    return Fraction(n1 if n0 * n1 < x else n0, 1 << s)

An alternative implementation of the original function could then be as follows:

def sqrt(x: Union[int, float, Fraction], precision: int = 53) -> Fraction:
    a, b = x.as_integer_ratio()
    precision += 3
    return sqrt_of_int(a, precision) / sqrt_of_int(b, precision)

This version has the advantage that sqrt(a/b)*sqrt(b/a) == 1, which the original does not guarantee. On the other hand, I'm not sure if the rounding guarantees hold anymore. But the error bounds are still guaranteed after increasing the precision by 3.

Overuse answered 15/10, 2023 at 3:13 Comment(0)
C
0

Newton's Method (simple code)

This is the method suggested in Think Python, 2nd edition, page 67, and doesn't need any library. Newton's method takes a number a and returns its square root as follows:

y = (x + a/x) / 2

where x is an arbitrary estimation and y is a better estimation of a.

def sqr_root(a): 
    if a <= 0:  # check whether *a* is a positive number.
        print('Non-positive entry')
    else:
        epsilon = 0.0000001  # an arbitrary small number, for faster convergence. 
        x = a/2  # an estimation of square root. 
        while True:
            y = (x + a/x) / 2 
            if abs(y-x) < epsilon:  # check the difference between x and y to be small enough
                break
            x = y
    print(x)

P.S: the Newton's Method formula in Latex

Counteroffensive answered 28/12, 2023 at 21:42 Comment(2)
This is a good idea, but the implementation is frankly terrible. A) The function needs to return a value, i.e. replace print(x) with return x. B) In case of an error, a function shouldn't just print and continue, it should raise an exception, i.e. raise ValueError('Non-positive entry'). C) Why doesn't it handle 0? You could maybe just add a special case, if a == 0: return 0. D) Least importantly, the precision is not great, e.g. sqr_root(1) = 1.0000000464611474. If I do epsilon = 1e-15 instead, I get 1.0 exactly and in only 7 loops, but I'm not sure that's a universal solution.Apulia
In that book, return is covered in section "6.1 Return values" and raise in section "11.4 Reverse lookup".Apulia
E
-5

find square-root of a number

while True:
    num = int(input("Enter a number:\n>>"))
    for i in range(2, num):
        if num % i == 0:
            if i*i == num:
                print("Square root of", num, "==>", i)
                break
    else:
        kd = (num**0.5)  # (num**(1/2))
        print("Square root of", num, "==>", kd)

OUTPUT:-

Enter a number: 24
Square root of 24 ==> 4.898979485566356
Enter a number: 36
Square root of 36 ==> 6
Enter a number: 49
Square root of 49 ==> 7

✔ Output 💡 CLICK BELOW & SEE ✔

Output1

Output2

Eastward answered 7/5, 2022 at 17:21 Comment(8)
Why do a linear search for an integer solution when you could simply do math.isqrt(num) ** 2 == num? Or even simpler, kd = math.sqrt(num); kd.is_integer()?Apulia
This produces no output for num = 2 or belowApulia
That else block triggers too often. For example, with num = 5, it prints Square root of 5 ==> 2.23606797749979 three timesApulia
if part is for the number which have integer square root and the else part is for float number square rootEastward
See the output i have tryed to use for loop and prime number logic.Eastward
That output doesn't match that code. Try it.Apulia
@Apulia it matches the output ... :)Eastward
✔👆 see the image bcz if any mistake in alignment then error will occur most offent ..... :)👆✔Eastward

© 2022 - 2024 — McMap. All rights reserved.